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Eo  = 1 keV,  n^=  3.34  x 10  ^ ///cc  and 

I = 1.44  X 10  ///sec  (liquid  density,  borderline 

chemistry) 228 

Figure  6.15.  Time  sequence  of  radial  profiles  for  liquid  density, 
borderline  chemistry  case  for  (a)  H,  (b)  OH,  (c)  H«, 

(d)  0,  (e)  HO2,  (f)  O2  and  (g)  H2O2 7 229 


Figure  6.16.  Time  course  at  r = 0 for  T = 300  K, 

~ ^ keV,  n2=  3.34  x 10  ///cc  and 

1 = 1.44  X 10  ///sec  (liquid  density,  fast 

chemistry) 233 

Figure  6.17.  Time  sequence  of  radial  profiles  for  liquid  density, 
fast  chemistry  case  for  (a)  H,  (b)  OH,  (c)  H,,  (d)  0, 

(e)  HO2,  (f)  O2  and  (g)  H2O2 234 

Figure  6.18.  Snapshot  of  near  steady  state  profiles  for  all  species 
for  low  density,  slow  chemistry  case,  t = 1 sec.  steady 
state  pure  diffusion  limit  overlaid  for  only— diffusing 
species 238 


Figure  6.19.  Snapshot  of  near  steady  state  profiles  for  all  species 
for  low  density,  borderline  chemistry  case, 
t = 1 sec 
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Figure  6.20.  Snapshot  of  near  steady  state  profiles  for  all, species 
for  I atra,  slow  chemistry  case,  t = 3.5  x 10~^  sec.  . 

Figure  6.21.  Snapshot  of  near  steady  state  profiles  for  all  species 
for  1 atm,  borderline  chemistry  case, 

t = 3.5  X 10  sec 241 

Figure  6.22.  Snapshot  of  near  steady  state  profiles  for  all, species 

for  1 atra,  fast  chemistry  case,  t = 3.5  x lo"^  sec 242 

Figure  6.23.  Snapshot  of  near  steady  state  profiles  for  all  species 
for  liquid  density,  slow  chemistry  case, 

t = 1.8  X 10  sec 243 

Figure  6.24.  Snapshot  of  near  steady  state  profiles  for  all  species 
^or  liquid  density,  borderline  chemistry  case, 
t = 1.8  X 10  sec 244 

Figure  6.25.  Snapshot  of  near  steady  state  profiles  for  all  species 
for  liquid  density,  fast  chemistry  case, 

t = 1.8  X 10  sec 245 

Figure  7.1.  Results  of  DIFKT  for  the  low  density,  slow  chemistry 
case  (unpublished  work  of  Rio,  Dayashankar, 

Schlppnick,  and  Green)  for  (a)  neutral  species, 

and  (b)  ionic  species  and  thermal  electron 258 

Figure  7.2.  Results  of  DIFKT  for  the  low  density,  borderline 

chemistry,  case  (unpublished  work  of  Rio,  Dayashankar 

Schlppnick,  and  Green)  for  (a)  neutral  species, 

and  (b)  ionic  species  and  thermal  electron 259 

Figure  A.  1.  Pure  diffusion  time  course  at  ? = 0 for  n = 0, 
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Diffusioa  kinetic  results  are  calculated  via  a modification  of  the 
BOAT  code  of  Dash  and  Pergament,  which  was  designed  to  calculate 
turbulent  jet  entrainment  with  chemistry  in  rocket  nozzle  exhaust 
plumes.  These  results  are  explained  in  terms  of  diffusion  modification 
of  the  corresponding  pure  kinetic  system.  Results  fall  into  two  regimes 
(early  and  late  chemistry)  depending  on  the  ordering  of  kinetic  onset 
time  and  diffusion  convergence  time.  Dependence  of  trends  on  injection 
rate  and  medium  density  are  identified.  Results  are  critically  compared 
with  steady  state  results  which  incorporate  the  same  inputs. 
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CHAPTER  I 
INTRODUCTION 

Background  and  Motivation 

The  irradiation  of  gaseous  and  liquid  media  by  either  photons, 
protons  or  other  positive  ions  ultimately  involves  the  generation  of 
electrons  and  their  subsequent  interaction  with  the  medium.  This 
interaction  takes  place  within  a time  interval  referred  to  as  the 
'physical'  phase,  involving  excitation  and  ionization  within  10“^^  sec, 
and  dissociation  of  molecules  within  10“^^  sec  (Hart  and  Platzmann, 
Magee,  1963).  In  the  case  of  liquid  water,  ion  molecule 
reactions  are  pushed  back  into  this  time  phase  (<  10“^®  sec)  by  the 
factor  of  concentration.  In  the  gaseous  phase,  this  is  not  the  case. 
This  stage  is  followed  by  an  intermediate  stage  known  as  the  'physico- 
chemical phase,  which  primarily  marks  the  end  time  of  dissociation  of 
molecules  and  the  beginning  of  ion  molecule  reactions.  Diffusion  of  the 
electron,  l.e.  the  thermalized  electron,  begins  at  this  time. 

In  the  liquid  case  the  'chemical'  phase  starts  around  -10“® 

sec,  at  which  time  neutral  reactions  begin  and  molecular  diffusion 
begins  to  have  an  effect.  The  events  occurring  in  this  phase  have  been 
treated  by  various  researchers  in  terms  of  an  assortment  of  different 
diffusion  kinetic  models.  Most  noteworthy  has  been  the  work  of 
Kupperman  (see  Kupperman  and  Bedford,  1962a,  1962b;  Kupperman,  1974), 
Magee  (1957)  and  Magee  and  Chatterjee  (1978,  1980). 

In  order  to  model  the  'chemical'  phase  behavior,  however,  physical 
Inputs  deriving  from  the  earlier  'physical'  phase  are  required.  This 
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connection  between  the  processes  and  events  taking  place  in  these  two 
phases  is  generally  known  as  the  interface  between  radiation  chemistry 
and  radiation  physics.  The  identification  of  the  relevant  early  time 
events  and  the  accurate  calculation  of  their  cross  sections  have  been 
the  object  of  intense  effort — and  the  subject  of  vigorous  debate — over 
the  past  several  years  among  those  Involved  in  this  effort  (for  example, 
see  Douthat  (1983),  Inokuti,  (1983)). 

A significant  contribution  to  this  effort  has  been  made  by  Green 
and  coworkers  at  the  University  of  Florida.  With  particular  regard  to 
the  problem  of  electron  degradation  in  water  vapor  are  the  works  of 
Green  et  al.  (1971),  Olivero  et  al.  (1972)  and  Kutcher  and  Green 
(1974a),  the  latter  work  extending  to  the  case  of  liquid  water.  These 
inelastic  cross  sections,  along  with  fits  to  elastic  cross  section  data 
by  Porter  and  Jump  (1978)  have  been  compiled  and  updated  by  Zaider, 
Brenner  and  Wilson  (1982)  at  Los  Alamos  National  Labs.  On  the  basis  of 
these  cross  sections  it  is  now  possible  to  calculate  'chemical'  phase 
events  on  the  basis  of  "first  principles,"  rather  than,  as  had 
previously  been  the  case,  by  assuming  resonable  spatial  distributions 
and  working  backwards  from  known  experimental  g-values  (yields  per  100 
eV).  In  this  regard,  two  widely  divergent  approaches  have  been  taken  in 
the  calculation  of  these  later  time  events.  The  Oak  Ridge  group  (see 
Hamm_e^al^. , 1978;  Hamm  et  al. , 1983;  Ritchie  ^ _al^. , 1978;  Turner  et 
_a^. , 1980;  Wright  et  al.,  1983a  and  b)  have  constructed  a Monte  Carlo 
program  which  follows  the  electron  into  the  chemical  stage  and  actually 
calculates  chemical  reactions  stochastically.  Such  an  approach  can 
probably  be  considered  to  be  state  of  the  art  in  simulation. 
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^ approach  has  been  taken  by  Green  and  coworkers.  This 

approach  consists  of  following  the  electron  in  its  degradation  and 
multiple  scattering  during  the  physical  phase,  at  which  point  a 
statistical  profile  of  Inelastic  events  is  taken  from  which  a smoothed 
mean  spatial  and  spectral  density  of  inelastic  events  is  obtained.  This 
spatial  spectral  density  is  called  the  spatial  yield  spectrum,  (Jackman 
and  Green,  1979),  generalizing  on  the  concept  of  the  yield  spectrum 
(Green,  Garvey  and  Jackman,  1977),  which  is  a modification  of  the 
concept  of  the  degradation  spectrum  of  Spencer  and  Fano  (1954).  It  is  a 
nearly  universally  Invariant  quantity  for  all  gaseous  media,  possesses 
advantageous  mixing  and  scaling  properties,  and  from  it  the  spatial 
yields  of  states,  and  from  them  the  spatial  yields  of  the  various 
iiiitial  deposition’  species  follow  in  a simple  manner. 

Why  the  Continuous  Electron  Beam? 

This  thesis  takes  as  its  starting  point  the  work  of  the  previous 
thesis  by  Rio  (1983)  which  calculates  spatial  yields  for  water  vapor  via 
the  Monte  Carlo  technique  alluded  to,  and  inputs  them  into  a time 
dependent  diffusion  kinetic  calculation  for  continuous,  steady  electron 
beam  turned  on  at  t = 0. 

Rio,  in  his  thesis  (1983),  did  attempt  to  calculate  steady  state 
diffusion  kinetics  of  this  system.  However,  at  this  point  in  time, 
these  efforts  remain  incomplete — close  to  final  answer  but  beset  by 
questions  of  numerical  technique  and  reproducibility.  Results  of  this 
study  will  be  of  use  in  the  resolution  of  this  problem. 

In  previous  and  other  current  radiological  work  a short  pulse 
electron  beam  rather  than  a continuous,  steady  beam  is  considered. 
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Kupperman  and  Bedford  (1962a,  b),  assuming  low  injection  rates  (also 
known  as_dose  rates),  and  therefore  isolated  tracks,  calculated  for  the 
case  of  an  initial  pulse,  the  effects  of  which  die  out  after  its  passage 
through  that  spatial  element  of  the  medium.  Magee  and  Chatterjee  (1978), 
implementing  the  'prescribed  diffusion'  model,  account  for  dose  rate 
dependence,  that  is  to  say,  nonlinear  effects  due  to  the  chemistry. 

They  also  calculate  for  an  initial  pulse.  The  present,  ongoing  Monte 
Carlo  work  of  the  Oak  Ridge  group  (see  previous  section)  also  calculates 
for  an  initial  pulse. 

In  the  microscopic  realm,  the  calculations  carried  out  here  for  a 
continuous  beam,  may  not  be  the  most  relevant  configuration.  However, 
for  high  enough  dose  rate  so  that  track  overlap  is  significant,  this 
macroscopic  approach,  averaging  out  track  structure,  should  be  valid. 

In  the  microscopic  realm,  it  may  also  apply. 

Outline  of  Thesis 

In  Chapter  II  we  begin  by  identifying  the  spatial  yields.  As  is 
noted  there,  it  was  decided  to  redo  the  work  of  fitting  the  spatial 
yield  spectra.  It  is  found  that  the  spread  in  the  spatial  fall  off  rate 
over  spectral  energy  of  the  electron  (involved  in  the  event  counted) 
can  be  adequately  represented  by  a span  between  its  two  extreme  spectral 
limits.  The  g— values  for  various  states,  and,  summing  over  branching 
ratios  devised  by  Rio,  for  the  various  primary  species,  can  then  be 
divided  into  a low  and  high  spectral  component.  It  is  found  that  this 
high  spectral  component  contributes  proportionately  less  and  less  at 
higher  incident  energy,  and  that,  as  a result,  there  is  a decrease  in 
charge  separation  at  higher  incident  energy. 
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The  exponential  fail  off  seen  by  Rio  In  the  data  out  to  two  to 
three  decades  is  confirmed.  The  data  are  presented  graphically  so  that 
there  can  be  no  mistaking  what  their  trends  are — to  within  the  noise 
limit  of  the  Monte  Carlo  data.  It  is  found  that  for  Incident  energy 
around  or  below  1 keV,  the  contours  are  sufficiently  close  to  spherical 
that  the  spherical  approximation  is  adequate.  This  is  discussed  further 
in  Appendix  D,  where  steady  state  axisymmetrlc  diffusion  is  calculated. 

In  Chapter  III  we  review  some  basic  mathematical  results  for  pure 
kinetics  and  pure  diffusion.  Of  specific  interest:  dependence  of 

reaction  time  on  source  strength,  reaction  rate  coefficient,  and  medium 
density,  and  dependence  of  diffusion  convergence  time  on  medium 
density.  We  then  address  diffusion  kinetics  of  two  body  decay.  First 
we  calculate  the  simplest  case,  namely  that  of  one  species  on  Itself. 

We  follow  this  up  with  the  case  of  that  between  two  distinct  species. 

In  the  two  species  case,  the  specific  reaction,  H + OH  + M > H2O  + M, 
the  context  of  the  radiological  deposition  is  used. 

The  two  species  diffusion  kinetic  problem  for  the  particular 
reaction,  H + OH  + M H2O  + M,  is  considered  at  the  intermediate 
density  of  the  three  medium  densities  considered  in  Chapter  VI,  at  all 
three  injection  rates.  In  this  way  its  effect  can  be  isolated.  We  find 
that  the  resultant  trends  can  be  explained  in  terms  of  (1)  the  ordering 
of  reaction  time  and  diffusion  convergence  time,  (2)  majority  and 
minority  reactants  (per  source  term),  and  (3)  relative  diffusion 


coefficients. 
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At  this  point  the  time  dependent  diffusion  kinetic  code  written  for 
-the  purpose  of  this  thesis,  DIFKIN,  is  brought  into  play.  This  code  is 
discussed  in  complete  detail  in  Appendices  A and  B,  In  effect. 
Appendices  A and  B are  a user's  guide. 

In  Chapter  IV  the  ion  reactions  are  considered.  From  a simple 
argument  we  identify  the  concentrations  of  neutral  species  subsequent  to 
local  electron-proton  recombination.  Noting  the  decrease  in  change 
separation  at  higher  incident  energy,  it  follows  that  at  higher  incident 
energy,  the  approximation  of  global  recombination  becomes  valid.  The 
possibility  of  upgrading  the  code  to  take  into  account  charge  separation 
exists,  and  should  prove  to  be  a very  interesting  direction  in  which  to 
take  this  work. 

This  chapter  is  closed  out  with  a calculation  of  dependence  of  the 
time  scales  of  the  ion  reactions  on  medium  density  and  injection  rate, 
based  on  results  of  Chapter  III.  This  calculation  is  based  on  a simple 
but  accurate  simplification  of  the  ion  reaction  set  (10  reactions)  into 
two  reactions. 

In  Chapter  V we  consider  the  pure  kinetics  of  the  neutral  reaction 
set.  Trends  in  medium  density  and  injection  rate  are  identified. 
Principal  pathways  and  their  trends  in  medium  density  and  Injection  rate 
are  identified  and  displayed.  It  is  found  that  above  a minimum 
injection  rate  qualitative  dependence  on  Injection  rate  disappears.  The 
question  of  stiffness  of  the  system  of  kinetic  equations  is  addressed. 

The  kinetic  systems  considered  are  those  which  obtain  at  the  center 
of  the  plume  for  the  representative  incident  energy,  1 keV,  ignoring 
diffusion.  It  serves  the  purpose  of  a base  for  comparison  so  that  the 
diffusion  effect  upon  kinetics  can  be  Isolated. 
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This  finally  leads  us  to  Chapter  VI  where  we  address  the  full 
diffusion  kinetic  system.  We  first  consider  the  trends  in  medium 
density  and  injection  rate  at  the  chosen  representative  medium  range 
incident  energy,  1 keV.  Specifically,  three  densities  are  calculated 
for:  6x10^^  #/cc,  1.84x10^^  #/cc  (1  atm  at  400  K) , and  3.34x10^^  ///cc 

(liquid  density).  The  first  two  cases  are  calculated  for  T = 400  K,  the 
last  one  for  T = 300  K. 

This  choice  is  made  for  the  sake  of  physical  realism.  All  kinetic 
runs  in  Chapter  V but  one  are  also  made  at  T = 400  K.  A corresponding 
run  is  made  at  300  K for  the  liquid  density  case  to  test  for  temperature 
sensitivity,  and  to  provide  a corresponding  basis  for  comparison  between 
kinetic  and  diffusion  kinetic  results  at  this  density. 

At  each  of  these  densities  the  trend  in  injection  rate  is  observed 
by  calculating  for  three  values  (exception:  two,  in  the  case  of  the  low 
pressure  vapor)  corresponding  to  late  chemistry,  borderline  chemistry, 
and  early  chemistry. 

Finally,  in  Chapter  VII,  the  results  obtained  here  are  compared 
with  results  of  Rio,  Dayashankar,  Schippnick,  and  Green  (unpublished 
work)  for  the  steady  state  case.  Future  directions  for  this  work  are 
discussed.  Experimental  set-ups  in  which  these  results  may  be  tested 
are  discussed,  especially  in  the  context  of  current  radiation  chemistry 
research,  in  particular,  the  Febetron  work  of  Hanrahan  and  coworkers 
(this  campus). 


CHAPTER  II 
SPATIAL  YIELDS 

In  previous  work  in  radiation  chemistry  predicated  upon  a 
macroscopic  diffusion  kinetic  model  (Kupperman  and  Bedford,  1961a, 

1962b;  Magee  and  Chatterjee,  1978,  1980),  the  "initial"  spatial  yields 
have  in  general  been  estimated  on  the  basis  of  known  experimental  g— 
values  (yields  per  100  eV),  assuming  a reasonable  spatial  dependence. 
Calculation  of  these  initial  spatial  yields  from  first  principles 
requires  (1)  a comprehensive  compilation  of  collision  cross  sections, 
and  (2)  solution  of  the  transport  equation  to  which  these  cross  sections 
are  inputs. 

A significant  contribution  to  this  first  task  has  been  made  over 
the  past  several  years  by  the  research  group  led  by  Alex  Green  at  the 
University  of  Florida.  Of  particular  interest  to  the  problem  of 
electron  degradation  in  water  vapor  are  the  work  of  Green  et  al.  (1971), 
Olivero  nt  al . (1972)  and  Kutcher  and  Green  (1976a),  the  latter  work 
extending  to  the  case  of  liquid  water.  These  cross  sections,  along  with 
fits  to  elastic  cross  section  data  by  Porter  and  Jump  (1978)  have 
recently  been  compiled  and  updated  by  Zaider,  Brenner  and  Wilson  (1983) 
at  Los  Alamos  National  Labs. 

In  his  Ph.  D.  thesis  (1983)  Daniel  Rio  has  recently  calculated 
these  initial  spatial  yields  by  means  of  a Monte  Carlo  technique 
developed  in  this  group  (Jackman,  1978;  Kutcher  and  Green,  1976b).  In 
this  calculation  the  primary  electron  and  any  secondaries  which  it  may 
have  created  are  degraded  randomly  according  to  probabilities  determined 
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by  the  set  of  collision  cross  sections  Inputted  into  the  code.  The 
collisional  history  raw  data  are  not  used  to  directly  calculate  the 
spatial  yields.  Instead,  the  spatial  yield  spectra  approach  of  Green  et 
(Jackman  and  Green,  1979;  Green  and  Singhal,  1979)  is  employed. 

This  approach  is  based  upon  the  calculation  of  a quantity, 
U(E,p,z,E^),  called  the  spatial  yield  spectrum,  or  four  dimensional 
yield  spectrum,  which  is  a function  of  incident  energy,  E^,  spectral 

energy,  E,  and  cylindrical  spatial  coordinates,  p and  z.  It  is  defined 
by 


U(E,p,z,E^)  = N(E,p,z,E^)/Ax  (2.1) 

where  N(E,p,z,E^)  is  the  number  of  inelastic  events  taking  place  per 
unit  time,  normalized  to  injection  rate,  with  spectral  energy  in  the 
interval  AE  centered  at  E,  in  the  volume  element  Ax  = 2Trpdpdz  centered 
at  (p,z).  Like  the  two  dimensional  yield  spectrum  (Green,  Garvey  and 
Jackman,  1977)  which  it  reduces  to, 

U(E,E  ) = U (E,E  )0(E  -E  -E)  + 6(E  -E)  (2.2) 

o a o om 

the  four  dimensional  yield  spectrum  is  divided  into  a continuous 
spectrum  extending  from  E^  up  to  E^  - E^,  U^(E, p, z,E^) , and  the 

contribution  of  the  undegraded  electrons,  called  the  source  term, 
D(p,z,E^)  (Green  and  Singhal,  1979): 


U(E,p,z,E  ) = U (E,p,z,E  )9(E  -E  -E)  + D(p,z,E  )6(E  -E) 
^ a o om  o o 


(2.3) 
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where  is  the  threshold  of  inelastic  events,  0(E)  is  the  Heaviside 
step  function  and  6(E)  is  the  Dirac  delta  function. 

Defining  p^(E)  as  the  probability  of  exciting  the  i-th  inelastic 
state,  given  an  inelastic  event,  the  spatial  yield  for  this  state, 
J^(p,z,E^),  is  given  by 


E 

J^(p,z,E^)  = / ° p^(E’)U(E',p,z,E^)dE'  (2.4) 

E . 

1 

where  E^  is  the  i-th  threshold  energy.  The  probability  p^(E),  is  given 
by  spatial  yield  spectrum  is  therefore  seen  to 

unite  into  one  function,  in  effect  a source  function,  the  common,  shared 
spectral  and  spatial  information  of  all  inelastic  events. 

Finally,  summing  over  the  appropriate  states  weighted  by  the 
appropriate  branching  ratios,  one  obtains  the  spatial  yields  of  the 
"primary  yield"  species  produced  by  the  degradation  of  the  incident 
electrons  (primary  electrons)  and  secondaries  and  tertiaries  (higher 
generations  are  negligible).  For  water  vapor  these  species  are  H,  OH, 

H2,  0,  H20''',  H"^,  OH"^  and  the  thermal  electron.  O'*"  production  (less  than 
1%)  is  negligible. 


The  Deposition  Spectrum 

Singhal  et  al.  (1983),  in  their  work  on  electron  energy  degradation 
in  neon,  define  a quantity  which  they  call  the  deposition  spectrum.  As 
the  electron  falls  below  threshold  it  is  counted  in  the  spectral  and 
spatial  bins  in  which  it  falls,  and  its  subsequent  elastic  degradation 
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is  not  followed.  Since  the  threshold  of  inelastic  events  in  neon  is 

fairly  high  (16.7  eV),  there  is  room  in  this  interval  for  spectral 
structure. 

The  low  inelastic  threshold  of  water  (threshold  of  rotational  and 
vibrational  states  is  .12  eV  precludes  this,  although,  if  one  wanted  to, 
one  could  keep  track  of  the  spectral  structure  below  the  threshold  of 
electronic  excitations  (7.6  eV).  Thus,  in  the  case  of  water  vapor,  it 
suffices  to  count  in  one  bin  below  threshold,  which  is  what  Rio  (1983) 
has  done. 

Thus,  apart  from  a normalization  constant,  this  deposition  spectrum 
for  the  case  of  water  vapor,  somewhat  in  the  manner  of  the  source  term 
at  the  other  end  of  the  spectral  scale,  is  merely  a scaled  spatial 
density  its  spectral  aspect  reduces  to  a mere  formality.  In  fact  Rio 
defines  the  deposition  spectrum  for  H2O,  DS(p,z,E^),  as  the  spatial 
density  of  subthreshold  electrons.  It  therefore  integrates  out  to 
N(Eq),  the  electron  yield  due  to  ionization  in  addition  to  the  initial 
or  Injected  electrons  (primaries). 

The  motivation  for  keeping  separate  count  of  these  events  was  to 
identify  the  spatial  yield  of  the  thermal  electrons,  i.e.,  those 
electrons  which  have  fallen  below  threshold  (at  300  K,  the  thermal 
energy  of  the  electron  is  about  .03  eV).  However,  as  we  shall  see 
below,  this  was  not  really  necessary.  At  these  low  energies,  the 
electron  has  essentially  stopped  moving,  and  its  spatial  distribution, 
to  within  a normalization  factor,  has  converged.  Thus,  to  within  a 
normalization  factor,  the  deposition  spectrum  in  the  case  of  water  vapor 
is  identical  with  the  low  spectral  end  of  the  spatial  yield  spectrum. 
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Analytic  Representation 


Green  and  Singhal  (1979)  represented  the  two  dimensional  yield 
spectra  calculated  by  Green,  Garvey  and  Jackman  (1977)  and  Green, 
Jackman  and  Garvey  (1977)  by  the  following  common  form: 


Eo/1000.  They  found  that  they  could  simultaneously  fit  yield  spectra 
for  a variety  of  gases  with  12  set  to  .585  and  L = 1 eV. 

The  first  term,  which  is  in  general  a small  constant  value,  serves 
to  represent  the  small  limiting  value  at  high  spectral  energy  which  is 
essentially  independent  of  incident  energy.  The  second  and  third  term 
together  represent  the  increasing  value  of  the  yield  spectra  with 
decreasing  spectral  energy. 

In  subsequent  work  Singhal  and  Green  (1981)  and  Singhal  et  al. 
(1983)  fit  the  spatial  yield  spectra  of  0 and  Ne  with  the  following 
generalization  of  the  two  dimensional  form  devised  by  Green  and  Singhal 
(1979),  called  the  microplume  model: 


U(E,E^)  = + C^X  + 


(2.5) 


where  X X(E,E^)  - Ej^/(E  + L)  with  E^  and  E measured  in  eV,  and  Ej^ 


U (E,p,z,E  ) = Z C.X 
3 o . „ 1 

i=0  R 


3 . A.(E  ) 

_ „ X o 


G^(p,z,E^) 


(2.6) 


A (E  ) 
s o 


(2.7) 


where 


Gi(p,z,Eo)  = exp[-  (■ 


a.  r 
1 


1+6 . z 


1 


(2.8) 
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and  the  coefficients  A^(E^)  are  normalization  factors  for  the  spatial 

integrals  of  the  G.(p,z,E^).  The  quantity  R(E^)  is  the  range  in  units 
2 

of  gm/cm  . The  spatial  yield  spectrum  is  in  units  of  #/eV/(gm/cm^)^. 
The  spatial  coordinates  are  in  units  of  R(E^).  The  quantities  G.  (and 
therefore  A^)  were  defined  in  a piecewise  fashion  for  the  forward  and 
backward  plumes,  i.e.  for  z > 0 and  z < 0.  Although  this  was  done  for 
the  sake  of  computational  convenience,  and  allowed  for  the  easy 
calculation  of  forward  and  backward  spatial  integrals,  it  results  in  an 
unphysical  discontinuity  in  the  spatial  contours  at  the  z = 0 plane. 

Rio  (1983)  observed  that  the  data  exhibited  essentially  exponential 
falloff  in  all  directions  outward  from  the  forward  peak,  and  elliptical 
contours  centered  about  the  peak.  He  devised  the  following  continuous 
form  to  represent  these  features: 


G^(p,z,E^)  - exp(-g^di/p^  + 6^(z-z^)^)  (2.9) 

where  d is  medium  density  in  gm/cc  and  p and  z are  in  units  of  cm.  As 
It  stands,  this  form  contains  one  mistake,  but  a reparable  one:  z 

o ’ 

which  is  presented  as  independent  of  d,  should  be  set  inversely 
proportional  to  d.  The  Rio  form  is  therefore  correctly  stated  by 

G^(p,z,E^)  = exp(-0^/p'2  + (2.10) 

where  p - pd  and  z’  = zd  and  z^  is  in  units  of  gm/cm^.  Knowing  the 

densities  at  which  the  various  runs  were  made,  the  correction  follows 
Immediately. 
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It  was  the  original  Intention  of  this  thesis  to  use  Rio's  fits  as 
inputs.  However,  other  mistakes  have  been  made  in  this  work,  and  for 
this  reason  the  author  decided  to  refit  these  data.  At  the  instigation 
of  Professor  Green,  a group  effort  was  then  undertaken  to  produce  data 
at  higher  energies  and  one  lower  energy  in  order  to  exhibit  trends  in 
incident  energy  more  clearly.  The  availability  of  these  extensions  to 
the  data  base  was  extremely  helpful.  The  author  is  especially  grateful 
to  Dr.  Dayashankar  and  Mr.  Jerome  Schwartz  for  their  work  in  the 
generation  of  these  data  sets. 

The  Yield  Spectrum 

The  two  dimensional  spectra,  as  Green,  Garvey  and  Jackman  (1977) 
observed,  possess  a fairly  characteristic  pattern  as  a family  of  curves 
when  plotted  against  E,  namely,  a fall  off  from  a flat  initial  level, 
scaled  very  nearly  by  E^,  which  falls  roughly  as  between  the  inverse 
first  and  second  power  in  E,  to  a constant  level  Independent  of  E at  the 
upper  spectral  end.  The  flat  initial  level  reflects  the  drop  off  of  the 
total  inelastic  cross  section  with  decreasing  energy  below  its  peak.  As 
previously  noted,  this  behavior  was  represented  by  Green  and  Singhal 
(1979)  by  the  form  given  in  Eq.  (2.1). 

If  the  yield  spectra  are  calculated  down  below  the  threshold  of 
electronic  excitations  into  the  vibrational  and  rotational  spectrum,  a 
rather  predictable  bookkeeping  result  follows.  Because  these  processes 
can  degrade  the  electron’s  energy  only  in  small  amounts,  the  frequency 
of  Inelastic  events  increases  as  the  electron  degrades  through  this 
region,  and,  correspondingly,  the  yield  spectra  start  to  increase 
steeply  with  decreasing  E in  this  region. 
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However,  the  details  of  the  yield  spectrum  In  this  low  range  are  of 
no  significance  to  the  main  purpose  of  the  yield  spectrum,  which  is  the 
calculation  of  spatial  species  yields.  These  low  lying  vibrational  and 
rotational  states  do  not  contribute.  We  therefore  limit  our  attention 
to  spectral  energies  above  .5  eV. 

The  original  form  of  Green  and  Singhal  is  found  to  well  represent 
the  two  dimensional  yield  spectrum  at  all  interesting  spectral 
energies.  In  fact,  it  can  even  be  simplified  without  loss  of 
accuracy  if  the  form  X = Ej^“/(E  + L)  is  replaced  by  the  form 

Q ~~2  2 

X = //  E + L , and  n is  fixed  at  .5  rather  than  .585,  and  L is 

to  vary.  can  then  be  set  to  zero,  and  we  obtain  = .0282, 

C2  = 902.,  and  L = 3.21  eV.  In  effect,  the  middle  plume  is  dispensed 
with.  The  two  dimensional  yield  spectrum  and  the  fit  to  it  are  graphed 

in  Fig.  2.1.  This  work  Involved  some  collaboration  with  Professor  Green 
and  Dr.  Dayashankar. 


The  Spatial  Yield  Spectrum 

Presentation  and  characterization  of  data.  In  considering  the 
spatial  yield  spectral  data,  since  the  two  dimensional  spectrum  has 
already  been  tied  down,  it  is  natural  to  normalize  it  out,  so  that  the 
remaining  spectral  dependence  stands  out  as  a separate  entity.  We 
therefore  define  normalized  spatial  yield  spectra  and  normalized 
deposition  spectra: 


U^(E,p,z,E^)  = U^(E,p,z,E^)/U^(E,E^) 


(2.11) 


(A3///)  (°a‘a:)n 
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Figure  2.1.  Two  dimensional  yield  spectra,  U(E,E  ) vs.  E. 
Data  and  fit.  ° 
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DS^(p,z,Eq)  = DS(p,z,E^)/N(E^).  (2.12) 

The  source  term  is  already  normalized  (see  Eq.  (2.7)). 

It  may  seem  strange  to  speak  of  "normalized  spectra”  in  that  what 
has  been  normalized  out  of  them  is  exactly  what  makes  them  spectra, 
i.e.,  namely,  that  they  should  reduce  to  spectra  upon  spatial 
integration.  What  they  are,  really,  is  spectrally  parametrized  families 
of  spatial  densities.  However,  the  term  "normalized  spectra",  even 
though  it  is  somewhat  anomalous,  is  convenient. 

In  order  to  accurately  characterize  these  data,  graphs  of  p and  z 
profiles  were  plotted  for  all  incident  energies  at  which  data  existed. 
Contour  plots  were  then  inferred  from  these  profiles  for  representative 
cases.  Representative  profiles  of  the  raw  data  at  E^  = 100  eV,  1 keV 
and  10  keV  are  shown  in  Figs.  2. 2-2. 4 for  the  normalized  deposition 
spectra  and  source  term,  and  in  Figs.  2. 5-2. 7 for  the  normalized  spatial 
yield  spectra.  Smoothed  contour  plots  are  also  shown  for  the  normalized 
spatial  yield  spectra  at  those  same  energies  in  Figs.  2.8—2.10. 

From  these  profiles  and  contour  plots  we  observe  the  following 
trend  in  E^:  for  low  E^  the  pattern  of  contours  is  concentric 

elliptical;  at  intermediate  E^,  the  pattern  becomes  singly  confocal; 
while  for  E^  > 10  keV  there  has  evolved  a highly  prolate  singly  confocal 
core  surrounded  by  a more  rounded  ’penumbra*. 

Turning  our  attention  next  to  the  100  eV  profiles,  we  note  that  the 
source  term  presents  a much  steeper  falloff  than  does  the  more  relaxed 
deposition  spectrum  (note:  the  source  term  is  only  noticeable  at  and 
below  500  eV,  no  data  exist  for  it  above  1 keV  for  the  number  of  primary 
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Figure  2.2.  Normalized  deposition  spectrum  ( ) and  source  term  ( ).  E = 100  eV. 
Data  and  fit.  ° 
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Figure  2.3.  Normalized  deposition  spectrum  E = 1 keV.  Data  and  fit 
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Figure  2.4.  Normalized  deposition  spectrum.  = lo  keV.  Data  and  fit 
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Figure  2.5.  Normalized  spatial  yield  spectrum.  E = 100  eV.  Data  and  fit 
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Figure  2.6.  Normalized  spatial  yield  spectrum.  E = 100  keV.  Data  and  fit. 
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Fipure  2.7.  Normalized  spatial  yield  spectrum.  'E'.  = 10  keV.  Data  and  fit. 
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Figure  2.8.  Smoothed  contours  of  rav7  data  of  the  normalized  spatial  yield  spectrum 
for  E 5 eVy  E = 10  eV,  and  E = 80  eV.  E = 100  eV. 
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Figure  2.9.  Smoothed  contours  of  raw  data  of  the  normalized  spatial  yield  spectrum 
for  E - 500  eV,  E = 800  eV,  and  E = 910  eV.  E =1  keV. 
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Figure  2.10.  Smoothed  contours  of  raw  data  of  the 
normalized  spatial  yield  spectrum  for 
E s 8 keV.  E =10  keV.  Pattern  is 
essentially  converged. 
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histories  implemented  in  the  generation  of  our  data  sets).  We  note  that 
its  peak  is  essentially  aligned  with  that  of  the  deposition  spectrum. 

Comparing  now  the  profiles  of  the  normalized  spatial  yield  spectra, 
it  becomes  clear  that  they  lie  in  between  the  limits  of  the  normalized 
deposition  spectra  and  the  source  term.  Contour  peaks  are  all  seen  to 
be  aligned.  Spectral  dependence  of  contour  shape  is  seen  to  be  none  at 
100  eV,  slight  at  1 keV  and  at  10  keV.  Fall  off  rate  decreases 
monotonically  from  a maximum  value  at  the  source  term,  and  the  contours, 
in  general,  are  seen  to  converge  to  the  contours  of  the  normalized 
deposition  spectrum  by  roughly  10  eV  for  = 100  eV,  by  roughly  500  eV 
for  Eo  =1  keV  and  by  roughly  8 keV  for  E^  = 10  keV. 

Theoretical  aside.  That  they  should  do  so  is  predicted  by 
continuous  slowing  down  theory.  The  stopping  power  in  the  Born  limit, 
dE/dx  cc  In  E/E  oc  (Bethe,  1932).  The  distance  traveled  up  to 

degradation  of  the  electron's  energy  down  to  E,  which  we  may  denote  as 
r(E,E^),  is  given  by 


r(E,E  ) 
o 


E 

f ° dE 
■'E  -dE/dx 


(2.13) 


From  this  it  follows  that  the  ratio  r(E,E^)/r(0,E^)  is  given  by 

r(E,E^)/r(0,E^)  = 1 - (E/E^)^/^  (2.14) 


from  which  it  follows  that  dr/dE  falls  to  zero  as  E approaches  zero 
The  quantity  r(0,E^)  is  defined  as  the  range  or  effective  range  and 
symbolized,  R(E^). 
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The  above  derivation  is  valid  if  is  well  enough  above  the  Born 
limit.  For  incident  energies  in  the  low  end  of  the  Bom  limit  and  below 
(typically  < 200  eV)  the  correction  to  the  stopping  power  based  on 
corresponding  correction  to  cross  sections  it  is  based  on,  given  by 
Green  and  Barth  (1965)  and  further  developed  by  Stolarski  and  Green 
(1967)  must  be  taken  into  account. 

Specifically,  the  stopping  power  mimics  the  cross  sections  of  which 
it  is  a weighted  sum,  and  turns  around  and  drops  to  zero  as  the  electron 
energy  drops  to  threshold.  As  a result,  the  energy  at  which  the 
electron  has  converged  in  its  motion  is  a smaller  fraction  of  incident 
energy  at  lower  incident  energy,  until  at  the  very  lowest  incident 
energies  the  convergence  has  disappeared.  This  relative  falling  back  of 
the  range  convergence  energy  with  decreasing  Incident  energy  is  shown  in 
Fig.  2.11,  where  r(E,E^),  inferred  from  these  spatial  yield  spectral 
data  by  means  of  continuous  slowing  down  theory,  is  plotted.  This  point 
is  discussed  in  greater  detail  below. 

Analytic  representation.  Three  basic  features  of  the  contours  at  a 
given  spectral  energy,  E,  can  be  Identified:  fall  off  rate  (in  essence, 
range),  offset,  and  shape.  In  the  interest  of  simplicity,  the  more 
differentiated  shape  of  the  contours  at  higher  Incident  energy  will  not 
be  addressed,  and  only  an  average  fit  will  be  attempted  there.  The 
following  form  is  chosen  to  represent  the  evolution  of  shapes  just 
described,  closely  at  low  and  medium  incident  energies,  and  in  an 
average  way  at  high  incident  energy: 


e(z ' 


G(p,z,E^)  = exp[-B(/y^p'^  + (z' 


(2.15) 


arbitrary  ualta 
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Figure  2.11 

Figure  2.11.  r(E,E^)  vs.  E for  representative  sequence  of 

values  of  E^.  Inferred  from  3(E^)  via  CSDA. 


in 
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for  which  the  normalization  factor,  A^(E^),  is  given  by 

A(E^)  = (l-e^)^Y^6^/8Tr  (2.16) 


This  form  may  be  described  as  affine  singly  confocal.  At  low 
energy,  e is  kept  equal  to  zero,  with  y > 1,  yielding  concentric 
prolate  elliptical  contours.  In  the  transition  range  y decreases 
towards  1 as  e takes  over.  At  higher  energies  e continues  to 
increase,  but  y decreases  below  1 in  order  to  represent  in  an  average 
sense  the  rounded  penumbra. 

The  true  eccentricity  of  the  elliptical  contours  is  given  not  by 

2 2 

e,  but  by  = /I  - (l-e  )/-y  . On  the  basis  of  the  stochastic  data 

obtained  so  far,  there  is  no  definitive  indication  of  a trend  toward 

oblateness.  We,  therefore,  constrain  e to  fall  back  to  zero  as  E 

cr  o 

increases.  In  this  way  the  roundedness  of  the  penumbra  is  represented 
at  the  expense  of  the  prolateness  of  the  inner  core. 

The  results  of  fitting  the  parameters  g,  z^,  e and  y to  the 
deposition  spectra  energy  by  energy  are  presented  in  Figs.  2.12-15.  The 
following  forms  were  found  to  adequately  represent  these  trends,  taking 
into  account  the  above-mentioned  nonperturbative  constraint  on  y: 


8(Eo)  = 


1+bE 


q,  q, 

+ C.E  ^ 

o o 1 o 2 o 


P ln(E  /E  +1) 

e(E  ) = r —2 — § ]Q 

o^  Lp  ln(E  /E  +1)  + iJ 
o e 


(2.17) 

(2.18) 


(2.19) 
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Figure  2.14.  e(E^)  vs.  E^.  Values  obtained  from  fitting  spectral 

data  at  individual  E values,  and  regression. 


Figure  2.15.  y(E^)  vs.  E^.  Values  obtained  from  fitting  spectral 

data  at  individual  E^  values,  and  regression. 
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-E  /E 


Y(E„)  - [1  + <,(E^E^)e  » 

o a 


a + e 


(where  y = /I  - e"^(E  )) 
00  o 


(2.20) 


where  the  energy  dependence  fitting  parameters  are  given  in  Table  2.1. 
The  quality  of  fit  is  indicated  in  Figs.  2.2-4  where  the  graphs  of  the 
fit  have  been  overlaid  on  the  data. 

The  form  2(E^)  is  identical  to  the  form  of  the  inverse  range  in 
Singhal  and  Green  (1981).  We  can  compare  our  results  with  Singhal  and 
Green,  and  with  experimental  data  of  electron  penetration  of  arbitrary 
media  of  Katz  and  Penfold  (1952).  These  data  and  more  recent  data  have 
been  fitted  by  Green  (private  communication,  1984)  by  the  formula, 

"'V  ■ (2.21) 

k 

where  = E^/1000.  The  corresponding  inverse  range  is  plotted  along 
with  B(E^)  in  Fig.  2.12.  It  is  to  be  noted  that  the  accuracy  of  the 
form  of  Eq.  (2.21)  is  less  at  lower  energies.  The  agreement  in 
proportionality  of  the  two  quantities  is  very  good,  with  a factor  of 
proportionality  of  about  10.  For  exponential  fall  off  this  implies 
setting  the  range  at  the  distance  at  which  the  spatial  integral  has 
included  99.5%  of  the  total.  This  is  a little  further  out  than  the 
usual  definition  of  range  (see  Green,  1955)  but  about  right. 

The  fit  to  8(E^)  can  be  improved  a little  by  going  over  to  the 

form 


S(E^) 


a 

E bE  P 
o o 


(2.22) 
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6(E  ) and  6 
o o 

a = 5.46(6), 
a = 1.34(7), 
k = 4.56 

z'(E  ): 
o o 

conf ocal 
Cj  = 7.29(-9) 

C2  = 1.11(-12) 

= 0.543 
q£  = 1>88 

e(E  ) : 
o 

E = 251. 
e 

P = 48.2 
Q = 90.6 


TABLE  2.1 

Fitting  Parameters  for  Spatial  Yield  Spectra 


b = 2.65(-5),  p = 1.66 

b = 2.65(-5),  p = 1.75 


.197 


Eq.  (2.17) 
Eq.  (2.18) 


concentric 
Cj  = 4.05(-9) 
C2  = 2.41(-12) 
q^  = 0.662 
q2  ~ 1.88 


Y(E  ): 
o 

ot  = 0.751 

E^  = 185. 

E = 407. 
a 

a = 736. 
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where  p 1.75  and  e is  a small  value.  Specific  values  of  a,  b and 
e are  given  in  Table  2.1.  This  form  is  preferable  to  the  form  of  Eq. 
(2.17)  because  it  satisfies  a requirement  of  continuous  slowing  down 
theory  (it  should  be  noted  that  Eq.  2.21  also  satisfies  this 
requirement)  - namely,  that  1 0 as  E^  10  (threshold  idealized). 

^ further  requirement  arises  from  the  previously  noted  correction 
to  the  energy  loss  function  (-dE/dx)  below  the  Bom  limit  which  states 
—dE/dx  turns  around  and  fall  to  zero  as  E^  decreases.  Noting  that 


1 


R’(e) 


(2.23) 


(this  result  follows  from  Eq.  (2.13)  and  the  definition  of  R(*)  via 
the  Fundamental  Theorem  of  Calculus),  it  follows  that  R'(E^)  must  turn 
around  and  rise  toward  infinity  as  E^  decreases  below  the  Bora  limit. 

The  form  of  Eq.  (2.22)  does  satisfy  this  requirement,  and  for  this 
reason  it  is  superior  to  the  form  of  Eq.  (2.17)  for  extrapolation  to 
lower  incident  energies.  It  is  on  the  basis  of  this  form  and  Eq.  (2.23) 
that  the  empirically  inferred  values  of  r(E,E^)  plotted  in  Fig.  2.11 
were  obtained. 

The  E^  dependence  of  z^  is  shown  to  be  approximately  proportional 
to  the  range,  as  expected.  The  concentric  z^  starts  to  increase  away 
from  the  true  z^  above  300  eV  to  compensate,  as  one  would  expect. 

Since  the  primary  reason  for  interest  in  the  spatial  yield  spectra 
for  the  purpose  of  this  thesis  is  to  use  them  to  obtain  spatial  yields 
for  diffusion  kinetics,  we  have  also  fit  the  spatial  yield  spectra 
with  e fixed  at  0,  i.  e.  using  the  concentric  form,  in  order  to  make  it 
easy  to  spherlcalize  them  in  order  to  use  them  as  inputs  in  the 
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diffusion  kinetic  calculations  carried  out  in  this  work  which  are  based 
upon  spherical  geometry.  Such  fits  will  be  somewhat  poorer  for  higher 
incident  energies,  because  of  the  inaccuracy  introduced  in  the  location 
of  the  offset.  However,  because  of  the  averaging  and  spreading  effect 
of  diffusion,  this  is  not  a serious  detraction.  The  full  benefit  of  the 
greater  realism  of  the  confocal  form,  or  of  a yet  more  apt  form,  can 
find  application  as  input  to  a general  two  dimensional  diffusion  kinetic 
code,  such  as  the  steady  state  code,  DIFKT,  written  by  Rio  (Ph.D. 
thesis,  1983). 

It  next  remains  to  fit  the  source  term  and  the  normalized  spatial 
yield  spectra.  All  that  remains  then  is  to  identify  the  tighter  fall 
off  rate  for  the  source  term  and  the  spectral  dependence  of  the  fall  off 
rate  over  the  spectral  range  of  the  normalized  spatial  yield  spectra. 

The  simplest  approach  would  be  to  represent  the  spatial  yield  spectra  by 
one  plume  with  6(E^)  replaced  by  B(E^,E),  and,  correspondingly,  A(E^) 
replaced  by  A(E^,E).  However,  this  disallows  a simple  procedure  for  the 
evaluation  of  the  spatial  yields  via  Eq.  (2.4)  in  which  the  spectral 
energy  dependence  has  been  separated  out. 

In  order  to  maintain  this  advantage  devised  by  Green  and  Singhal 
(1979),  we  gate  between  spectral  limits  instead  outside  the 
exponential.  Thus  we  set 

U^(E,p,z,E^)  = f(E,E^)A^(E^)G^(p,z,E^)  + 

(2.24) 

(l-f(E,E^))A2(E^)G2(p,z,E^). 

The  form  of  f(E,E^)  implied  by  the  Green-Singhal  form  is 
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f(E,E^)  = 


? 9 

Cq(E  + L^) 


Cq  + C^(E^  + L^)  + C^E^/IOOO 


(2.25) 


Contiuous  slowing  down  theory  suggests  a gate  of  the  form 


f(E,E^)  = 1 - r(E,E  )/R(E  ) 


e"  + bEP 

E ^ + bEP  ' 
o 


(2.26) 


We  may  adjust  all  three  parameters  e,  b and  p to  optimize  the 

fit.  We  now  have  two  candidate  spectral  gate  forms. 

It  remains  also  to  identify  S^(E^),  and  6g(E^).  We  may 

at  once  set  62(E^)  = S(E^).  Correspondingly  we  may  set 

^o^^o^  sake  of  clarity,  the  notation  3(E  ) will 

o 

be  replaced  by  the  notation  3/E^).  It  remains,  therefore,  to 
determine  3^(E^)  and  f(E,E^).  (it  also  remains  to  be  determined 
whether  the  form  of  Eq.  (2.24)  or  the  intra-exponential  gate  form 


= A(E  ,E)G(p,z,E  ,E) 


(2.27) 


where 


3 = 3(E^,E) 


f(E,Eo)3s(Eo)  + (1  - f(E,E^))3^(E^) 


(2.28) 


more  faithfully  represents  the  data. 

Fitting  the  source  term  data  (spectral  gate  independent)  and  the 

spectral  data  for  3^(E^),  we  find  that  we  obtain  good  fit,  using  the 
empirical  form 
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e (E  ) = — 5^2 

" ° E " + kbE  P 

o o 


(2.29) 


with  k - 4.56.  Thus,  the  faster  fall  off  of  the  high  spectral  plume 
becomes  progressively  relatively  less  so  with  Increasing  E^. 

Fixing  k,  and  varying  the  adjustable  parameters  of  the  CSDA 
suggested  gate,  we  obtain  a sevenfold  decrease  in  with  the  new 

values,  e’  = .674,  b'  = 2.99  (-10),  and  p'  = 3.91.  This  value  of  the 

2 . 

X IS  essentially  identical  to  that  obtaining  for  the  Green-Singhal 
gate! 

This  work  suggests  that  there  must  be  something  very  right  about 
this  empirical  form!  This  comes  as  something  of  a surprise,  and  makes 
this  pursuit  more  interesting. 

An  idea  of  the  quality  of  fits  can  be  gotten  from  viewing  Figs. 
2.16-18  where  contours  of  fits  have  been  plotted  for  representative 
spectral  plumes  — unfortunately,  for  the  unadjusted  parameters  in  Eq. 
(2.29)  and  for  the  wrong  k,  (k  = 3.6),  and  comparing  with  the  smoothed 
contours  of  the  raw  data  in  Figs.  2.5-7. 

It  might  be  commented  that  the  ideal  manner  in  which  to  proceed 
with  this  work  is  to  display  the  spectral  plume  data  in  interval  contour 
plots,  and  so  that  a precise  idea  of  the  shape  of  the  spectral  gate  can 
be  obtained. 


Spatial  Yields 

Spatial  yields  for  the  various  excited  states  (which  either 
dissociate,  ionize,  or  dissociate  and  ionize)  and  ionization  states 
created  by  the  electron’s  energy  deposition  are  calculated  via  Eq. 
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Figure  2.16.  Contours  of  fit  to  normalized  spatial  yield  spectrum  for  E = 5 eV,  E = 10  eV 
and  E = 80  eV.  E = 100  eV.  ’ ’ 
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Figure  2.17.  Contours  of  fit  to  normalized  spatial  yield  spectrum  for  E = 500  eV,  E = 800  eV 
and  E = 910  eV.  E = 1 keV. 
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Figure  2,18.  Contours  of  fit  to  normalized  spatial  yield  spectrum  for  E = 4 keV,  E = 6 keV, 
and  E = 8 keV.  E = 10  keV. 
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Figure  2.18.  Contours  of  fit  to  nomalized 
spatial  yield  spectrum  for 
E <_  8 keV.  E = 10  keV. 

Pattern  is  essentially  converged. 
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(2.4).  Spatial  yields  for  the  various  species  are  then  found  by  summing 
over  the  appropriate  branching  ratios.  By  virtue  of  the  device  of 
separation  of  spectral  energy  dependence  from  spatial  dependence 
mentioned  at  the  end  of  the  last  section,  the  integral  for 
'^£(p>z,E^)  simplifies  to 


J. (p,z,E  ) = 
1 o 


where 


(2.30) 


E -E. 

“lj‘V  ■ \<E')CjXJ(E’,E__)dE',  j . 0.2  (2.31) 

i 

and  ct.^(E^)  = p.(E^). 

Matters  are  made  even  simpler  by  our  identification  of  the  spatial 
dependence  of  the  source  term  with  the  spatial  dependence  of  the  high 
spectral  energy  plume,  i.e.  by  setting  G^(p.z.E^)  = G^(p.z.E^),  thereby 

reducing  the  number  of  separate  spatial  dependence  terms  in  Eq.  (2.30) 
to  two. 


Component  g— values 

The  coefficients,  a_(E^),  define  associated  g-values  via  the 
relation,  = g^^E^/lOO.  For  the  sake  of  notation  we  may  define 

“io,s  ~ “io  ^ “is  the  sum  of  these  two  quantities  is  of 

interest.  We  may  also  define  the  quantity  a.  = a + « . it 

i 12  io,s* 
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the  associated  g-value,  g.  which  is  generally  referred  to  as  "the"  g- 
value,  and  which  shall  be  referred  to  here  as  the  total  g-value.  Ihe 
idea  of  defining  component  g-values  follows  naturally  from  the  Green- 
Sxnghal  concept  of  separate  spectral  plumes. 

In  his  Ph.  D.  thesis  (1983)  Rio  calculated  the  total  g-values  of 
the  various  states  as  functions  of  and  displayed  them  graphically. 
From  these  results  via  the  branching  ratios,  the  g-values  for  the 
various  species  follow  without  too  much  trouble. 


The  total  g-values  (g^)  and  high  spectral  component  g-values 
(Sio,s>  graphed  here  in  Fig.  2.19  for  the  eight  individual  species: 
e , H2O  , H , OH  , H,  OH,  H2,  and  0,  and  for  the  sum  of  ions.  Note  that 

the  high  spectral  component  dies  off  as  increases  above  500  eV.  This 

reflects  the  increasingly  earlier  spatial  convergence,  spectrally 
speaking,  of  the  spatial  yield  spectrum  with  increasing  E^,  noted  and 
discussed  in  the  previous  section.  Specifically,  the  spectral  gate, 
f(E,E^),  as  a factor  in  the  integrand  of  the  probability  integral,  acts 
like  a tunable  high  pass  filter.  As  a result  of  this  tuning  over  E^, 

ratio,  gio,s'^®i’  ^ ^ ^o  decreases,  and  decreases  to  a small 
steady  level,  about  0.1,  at  higher  E^,  as  can  be  seen  from  Fig.  2.18. 

At  this  point  it  is  important  that  we  recall  our  empirical  result 

from  chi  square  fitting  of  the  data  and  graphical  display  that  the  short 

range  of  the  high  spectral  plume  relaxes  as  E^  increases.  If  this  is  a 

true  result,  it  means  that  the  effect  of  the  high  spectral  plume  on  the 

relative  profiles  of  the  spatial  yields  is  further  decreased  at  higher 
E_. 
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Charge  Balance 

Another  matter  comes  to  light  in  these  graphs.  This  is  the  matter 
of  charge  balance.  The  electron  yield,  other  than  for  contribution  from 
ionizing  decay  of  excited  states,  is  calculated  independently  of  the 
corresponding  ion  yields,  namely  directly  from  counting  in  the  Monte 
Carlo  calculation,  and  therewith  the  normalization,  N(E^),  in  the 
deposition  spectrum. 

In  these  graphs,  for  every  species  except  the  electron,  an 
uncorrected  and  a corrected  value  is  shown.  The  uncorrected  value  is 
calculated  according  to  prescription.  The  corrected  value  (i.  e. , 
corrected  in  the  sense  of  maintaining  the  constraint  of  charge  balance) 

is  calculated  according  to  a simple  renormalization  procedure  described 
below. 


Upon  comparison  of  the  g-value  of  ion  pairs,  obtained  by  summing 
the  g-values  of  all  three  ion  species,  with  the  g-value  of  ion  pairs 
implied  from  the  g-value  for  the  electron  (subtracting  out  the 
primaries),  a discrepancy  is  noted  which  is  small  at  low  energies  but 
which  Increases  from  200  eV  and  becomes  significant.  This  is  seen  in 
Fig.  2.20  in  which  the  corrected  g-value  for  ion  pairs  is  made  to  agree 
with  the  g-value  for  ion  pairs  implied  by  the  g-value  for  electrons. 

The  source  of  the  problem  is  not  clear.  It  could  be  a result  of 
misalignment  of  the  total  inelastic  cross  section  with  the  sum  of  the 
individuals~or  there  may  be  a problem  with  the  Monte  Carlo  program-one 
of  Inaccuracy  in  assignment  of  the  two  dimensional  yield  spectrum.  This 
latter  possibility  can  be  discounted  since  this  part  of  the  code  was 


Figure  2.19. 


Total  g-value  and  high  spectral  plume  component  for 
(a)  the  thermalized  electron,  (b)  positive  ions, 


(c)  (d)  H'^,  (e) 


OH' 


(f)  H,  (g)  OH,  (h)  H2 
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Figure  2.19  (a)  thermal  electron 


Figure  2,19  (b)  positive  ions 
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Figure  2.19 


(c) 
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Figure  2.19  (d)  H 


Figure  2.19  (e)  OH 


100000 


50 


Figure  2.19  (f)  h 


Figure  2.19  (g)  OH 
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Figure  2.19  (i)  0 
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checked  out  by  Green,  Garvey  and  Jackman  (1977)  against  the  modified 

discrete  energy  bln  (MDEB)  method  of  Peterson  (1969)  and  mas  found  to  be 
in  good  agreement. 

However,  this  problem  need  not  deter  us.  Granted  that  if  something 
is  amiss  with  the  ions,  there  may  be  reason  to  suspect  the  rest, 
however,  the  discrepancy  is  not  great:  10%  at  1 keV,  16%  at  10  keV. 
Therefore,  a little  shifting  around  of  values  does  not  do  any  harm  even 
at  worst,  and  is  probably  better  than  that.  The  important  thing  sought 
IS  exact  charge  balance— as  a constraint.  The  following  corrective 
procedure  therefore  seems  reasonable:  renormalize  the  ionization  states 

by  the  ratio  N(E^)  - 1 divided  by  the  total  ion  yield  per  incident 
electron  as  is. 

Since  the  discrepancy  in  ev/ion  pair  is  not  very  great,  the 
corrective  procedure  outlined  above  would  seem  to  be  adequate  to  our 
purposes  here.  As  can  be  seen  from  Fig.  2.19,  the  effect  of  the 
procedure  on  the  three  major  species  (other  than  the  electron,  which 
remains  unaffected),  namely  H,  OH,  and  H2O+,  is  about  the  same  as  for 

ev/lou  pair  In  the  cane  of  H3O+  (the  .vajor  ionic  apeclea),  aomewhat  less 
for  H,  and  nil  for  OH. 


In  Fig.  2.20  are  graphed  the  (corrected)  g-values  of  the  post 
recombinations  yields  assuming  the  approximation  of  globally 
instantaneous  recombination  of  ion  pairs  (see  Chapter  IV),  l.e.. 
Ignoring  charge  separation. 


Figure  2.20.  Total  post-recombination  g-value  and  high  spectral  plume 
component  for  (a)  H and  OH,  and  (b)  H2  and  0 assuming 
global  instantaneity  (see  Chapter  IV). 
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Figure  2.20  (b) 
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Charge  Separation 

The  next  matter  that  comes  to  our  attention  is  charge  separation. 
Noting  the  difference  between  the  large  high  spectral  plume  component  of 
the  g-value  for  positive  ions  and  the  small  value  for  the  electron  (this 
small  value  is  a post-Monte  Carlo  result,  and  arises  from  the  decay  of 
certain  excited  states  of  H2O  (see  Rio,  1983)).  it  is  clear  that  charge 
separation  is  significant,  and  more  so  at  lower  incident  energy.  A good 
illustration  of  this  is  provided  by  Kupperman  (1974)  where 
representative  profiles  of  H30^^,  OH,  and  e^^  are  plotted  for 
typical'  (?  eV)  Incident  energy  for  liquid  H2O. 

The  reason  for  mentioning  charge  separation,  other  than  to  call 
attention  to  one  imporant  feature  of  the  g-values  and  component  g- 
values,  IS  to  emphasize  the  sensitivity  of  the  accuracy  of  their 
calculation  to  the  accuracy  of  representation  of  the  yield  spectra. 
Specifically  what  is  meant  here  is  that  the  value  of  the  high  spectral 
plume  g-value  (taking  the  Green-Singhal  approach)  is  highly  sensitive  to 
the  spectral  variation  of  f(E,E^).  It  would  be  advantageous,  as  a 
future  direction  for  work  in  this  effort,  to  extend  the  Monte  Carlo  data 
reduction  to  include  numerical  calculation  of  spatial  yields  along  with 
spatial  yield  spectra  so  that  definitive  comparison  of  direct 

calculation  with  calculation  via  the  intermediary  of  the  fitted  spatial 
yield  spectra  would  be  possible. 


VJ-value 

It  is  of  interest  to  represent  one  particular  g-value  analytically, 
namely  the  low  spectral  plume  component  for  electrons,  g 

e d 
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representing  the  contribution  of  the  deposition  spectrum.  (Here  the 
notation  d rather  than  2 is  used  to  refer  to  the  low  spectral  plume, 
since  It  is  for  the  electron,  and  therefore  arises  from  the  deposition 
spectrum.  This  is  done  for  the  sake  of  formal  consistency.)  The 
quantity  more  frequently  referred  to  in  this  context  is  the  W-value,  or 
ev/ion  pair.  From  its  definition,  it  follows  that 


g _ - (100/E  )a  = (100/E  )N(E  ). 

e d ° e"d  o o 

We  may  define  g,p  = (100/E^)(N(E^)  - 1).  the  g-value  of  ion-pairs  which 
discounts  the  primaries.  It  follows  that  W = 100/g^p  and,  therefore,  W 

= VWE„)  - 1). 

Noting  that  N(E^)  must  fall  to  1 as  E^  decreases  to  I = 12.6  eV, 
the  ionization  threshold,  and  observing  the  near  proportional  LET  limit, 
we  can  conveniently  represent  N(E^)  by 

r 1 + 0.0415  (E  E > I 

N(E  ) = J ° o 

kl,  E < I 
o 


(2.32) 


CHAPTER  III 

ANALYTIC  PRELIMINARIES 
Introduction 

The  purpose  of  this  chapter  Is  to  gain  some  Insight  Into  the 
kinetics  of  both  the  ion  and  the  neutral  reactions  and  diffusion 
kinetics  of  the  neutral  reactions  by  abstracting  to  simpler, 
representative  cases.  It  is  possible  to  treat  many  of  these  cases 
analytically— either  exactly  or  approximately.  Numerical  results  are 
also  made  use  of,  and  are  obtained  via  a program  written  for  this 
thesis,  called  DIFKIN  (adapted  from  the  BOAT  code  of  Dash  and  Pergament 
—(1978),  see  Appendices  A and  B),  which  calculates  time  dependent 
diffusion  kinetics  in  spherical  geometry. 

In  the  case  of  pure  kinetics  both  the  case  of  an  initial  pulse  and 
of  a continuous  pulse  turned  on  at  t=0  are  discussed.  This  is  for  the 
most  part  for  the  sake  of  background,  since  these  basic  results  are  well 
known  (for  example,  see  Leu  et^.,  1973).  The  treatment  of  two  body 
distinct  species  decay  for  continuous  yield  is,  however,  original,  to 
the  best  knowledge  of  the  author. 

In  order  to  clarify  the  role  of  diffusion  in  diffusion  kinetics, 
the  case  of  pure  diffusion  is  first  presented  for  the  sake  of  background 
(in  spherical  geometry).  The  case  of  time  dependent  diffusion  in 
axisymmetric  geometry  is  discussed  in  Appendix  C,  although  this  work  is 
as  yet  incomplete.  Results  for  steady  state  diffusion,  however,  are 
presented  and  discussed. 
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For  diffusion  kinetics  only  the  continuous  pulse  is  discussed.  We 
first  consider  t„o  body  decay  of  one  species.  The  basic  £ln,e  dependent 
trends  are  identified  and  discussed.  The  asymptotic  limits  of  the 
steady  state  solution  are  also  identified  and  discussed.  An  analytic 
spline  between  these  limits  is  identified  and  presented  in  Appendix  D. 

The  case  of  two  distinct  species  is  considered  next  in  a physically 

specific  context,  namely  the  recombination  reaction,  H + OH  + M x H2O  + 

M,  in  order  that  we  might  be  able  to  Isolate  its  effect  in  the 

corresponding  calculated  results  for  the  full  system  later  in  Chapter 
VI. 

The  general  case  of  diffusion  kinetics  of  v-body  decay  is  also 
considered  as  an  aside  of  theoretical  Interest  in  Appendix  E.  The 
concept  of  an  Inverse  length  associated  with  the  reaction  rate  and 
diffusion  coefficient  is  generalised,  and  the  trend  with  v is 
investigated. 


IX  Lilt;  L L C S 


Single  species  decay.  We  begin  with  a simple  result  - the  decay 
of  a single  species  by  reaction  with  itself.  We  may  first  consider  the 
general  v-body  case  and  then  specialise.  Denoting  the  concentration  by 
c,  the  "yield  rate"  or  source  term  by  y,  and  reaction  rate  by  k.  the 
equation  in  the  pure  kinetic  case  is 


dc/dt  = -vkc^  + y 


(3.1) 


which  for  continuous 


source  (y  ^ 0)  becomes 
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dy/ciT  = 1 - y' 


(3.2) 


upon  ...Icing  ihe  transfor^attonn  ■,  . “AiET^  e and  i . ya/"/7AS. 
is  easily  solved  In  the  cases  u - 1 and  v - 2 to  give  respectively 


= - f (1  - e-“') 


(3.3) 


and 


' ~ ^ tanh^Zky  t. 


(3.4) 


In  the  general  case  the  Inverse  solution  Is  easily  found  by  expansion  in 
geometric  series  and  Integrating: 


T = Z 


1 


vj+1 


j=0  vj  + 1 


(3.5) 


The  trend  of  the  solution  In  v Is  shown  by  the  graph  In  Fig.  3.1  In 
which  Y(T)  Is  plotted  versus  t for  a sequence  of  values  of  v.  The 
trend  Is  as  one  would  expect  Intuitively,  namely,  that  as  the  number  of 
molecules  (or.  In  general,  fractions  of  units)  which  must  come  together 
to  form  a reaction  Increases,  the  approach  to  saturation  becomes  more 
abrupt  In  time.  This  Is  seen  already  1„  the  cases  v - 1 and  v - 2. 

The  value  of  y at  reduced  time  r . 1 (at  which  the  asymptotes 
Intersect)  Is  1-1/e  In  the  case  v - 1 and  1 - 2/(e2(l  . l/e2))  in  the 
case  V . 2.  In  general  for  v > 1 a 'reaction  time'  or  'transition 
time'  may  be  defined  at  x . i.  yhe  transition  Is  fairly  abrupt 
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between  beginning  and  end.  For  0 < v < 1 onset  is  gradual  ritlle  end 
of  transition  is  abrupt  (in  logarithm  of  time). 

The  initial  pulse  case  can  be  solved  explicitly,  and  is  given  by 


r 


(l  + v(v-l)c  ’ 

o 


V ^ 1,  V > 0 


c = V 


(3.6) 


-kt 

^ c^e  , V = 1 


These  solutions  are  plotted  in  Fig.  3.2  as  c/c  vs.  c '^^kt  for  a 

o o 

sequence  of  values  of  v. 

The  curve  for  v = 1 is  a tangential  envelope  of  those  for 
V > 1.  The  trend  with  v for  v > 1 is  analogous  to  that  in  the 
case  y ^ 0,  only  instead  of  more  abrupt  here  the  trend  is  towards  more 
gradual.  This  is  because  here  the  system  is  simply  decaying,  while  the 
decay  is  offset  by  the  competing  effect  of  the  continuous  source  term. 
Thus,  at  first,  decay  is  insignificant,  only  the  effect  of  the  source 
term  is  seen.  Finally  because  of  the  built  up  concentration,  decay  has 

gained  strength,  and  the  balance  shifts.  This  time  can  be  called  a 
'critical'  time. 

In  the  pure  decay  (initial  pulse)  case,  an  interesting  result 
obtains  for  v < 1.  That  is,  that  the  decay  accelerates  until  an 
endpoint  time  is  reached  at  which  the  material  is  completely  consumed. 

The  effect  of  stoichiometric  factor  is  irrelevant  to  this,  only  the 
power  law  matters. 
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Two  body  two  specte.  drcay.  We  next  consider  the  case  of  coupled 

two  body  decay  between  two  distinct  reactants.  The  kinetic  equations 
are  given  by 


dc^/dt  = -kc^c^  + 


dc^/dt  = -kc^C2  ^2 


(3.7) 


het  us  first  consider  the  Initial  pulse  case.  Setting  y,  - yj  - 0,  we 

obta.n  Cj  - Cjq  + (c2  - C2q).  The  equations  thereby  decouple  to 

give 

dc  /dt  = k(c,^  - c )c  - kr  ^ 

1 10  20^  1 

(3.8) 

dc„/dt  = -k(c,„  - c )c  - kr  ^ 

2 10  20^  2 ^2 


The  solution  is  found  by  integration  and  is  given  by 


^1  = 


^lO^^'lO  " ^^20^ 


‘^lO  ~ ^20® 


"‘^^‘^10  - ^20^^ 


^2 


.^2^10  ~ ^20^ 


^10^ 


^^*^10  “ ^20^^ 


(3.9) 


- c 


20 


This  result  is  in  agreement  with  intuition— the 
off  to  the  initial  difference  while  the  minority 
destroyed.  The  asymptotic  limit  (Eq.  (3.6),  v = 


majority  species  falls 
species  is  completely 
2)  is  approached 
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“ '10  - '20  * » O'.  "O'O  generally,  as  (c^^,  - * 0 If 

we  replace  c by  c/2  to  reflect  counting.  In  the  opposite  limit, 

*'10  ■ '2O^''^'l0  '20^  * *’  '■'O  asymptotic  limit  (Eq.  (3.6), 

w - 1)  Is  approached  for  C2,  with  k In  this  equation  replaced  by  kcjg, 
where  this  k Is  from  Eq.  (3.9).  Thus,  Cj  Is  effectively  a constant 
background,  while  Cj  undergoes  effective  one  body  decay.  This  equation 
was  first  Investigated  by  Verhulst  (1838)  In  the  study  of  saturation  of 
populations  (see  Green  and  Burkl,  1974). 

We  next  consider  the  continuous  pulse  case.  Making  the 

transformation  a,  = (c,  + Crvl/?  ana  =,  - / \ /n 

1 12  ^2  ~ ^2^  ecj^^tfons 

transform  to 


daj/dt  = -k(a^2  _ ^^2^  ^ ^ 


da2/dt  = (y^  - y^)/2. 


(3.10) 


The  second  of  these  equations  immediately  integrates  to  give 


^2  = 


(3.11) 


which  upon  substitution  into  the  other  equation  yields 


da^/dt  = -k(a^"  - (y^  - y2)2t2/4)  + (y^  + ^ 


(3.12) 


We  make  the  transformation,  x = kt,  yielding 


da 


l/dx  = -a^2^  x2.^ 


2k 


(3.13) 
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This  is  a first  order  dlffereetlal  equation  of  Rlcattl  type  for  »hlch 

the  transformation,  a^  = u'/u,  transforms  the  equation  to  a second  order 
linear  differential  equation.  We  obtain 

u"  = (B^x^  + C)u  (3.14) 

where  B - (y^  - y2)/2k  and  c = (y^  + y2)/2k.  Other  than  for  the  sign 
of  C,  this  is  just  the  Schroedinger  equation  for  a simple  harmonic 
oscillator.  Following  the  development  in  Eisberg  (1961),  we  may 
transform  one  more  time,  setting  = c/B  = (y^  + ~ Y^)  and 

5 X = /(y^-y2)k/2  t to  obtain 

u"  = (5^  + a^)u. 

The  solution  is  found  by  making  the  substitution  u 
the  differential  equation  for  H, 


(3.15) 


5^/2 

= He"*  ' , yielding 


H"  + 2CH'  + (1  - a^)H  = 0 (s.jg; 

for  which  the  solution  is  most  conveniently  found  by  means  of  a power 
series.  The  recursion  relation  is  given  by 

2 

h = « - 1 - 2n  . 

n+2  (n+l)(n+2)  n’  ^ ^ 0*  (3.17) 

The  Initial  condition  a,  - 0 Implies  that  hj  - 0;  Is  only  required 
not  to  be  zero,  and  is  arbitrary  as  it  scales  out. 
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Solving  now  for  we  obtain 


- (du/dx)/u  = /B  (du/d?)/u 


= /b(5  + H'/H) 


(3.18) 


(yi  - Y2)t  . y,  - y 


+ /-_L_2  H'(g;g) 

2k  H(5;a) 


^ " ’'(^1  - yo)k/2  t and  = (y  + V )/('„  „ s 

^ ^ VYj  1-  y2;/(.y^  - y^)  as  noted 

above.  Noting  that  hj/h^  . (a^  - „/2,  it  follows  that 

®l  ~ (yj  + y2>c/2  for  t ~ 0.  From  the  fact  that  H’/H  * 0 as  { a . (to 
easily  see  this,  consider  the  case  where  H truncates  to  a polynomial)  It 
follows  that  aj  t (y  - y ,3  3 __  Therefore  c ~ y t 


2 ^2  r t 0 and  ~ (y^  - y^)t  and  ^ 0 as  t ->■  ®,  as  we 


and  c 

might  Intuitively  expect 

From  this  we  see  that  the  quantity, 


gU:«)  - H’(C;o)/[(c.^  - l)CH(5;a)] 


gates  the  solution  between  Its  asymptotic  limits  In  t.  The  essential 

Information  of  the  solution  Is  therefore  seen  to  be  contained  In  this 

function  in  its  most  meaningful  form.  In  order  to  get  a better  look  at 

g(C;a)  we  may  define  an  alternative  dimensionless  time,  r,  - t,  and 

a new  parameter,  a . (y^  . ^ ^ 

of  5 and  d by  a . (o^  - l)/2  and  n . /I/S  5,  The  purpose  of  this  new 
variable,  „,  is  to  provide  a fixed  point  of  reference  while  varying  y, 
while  keeping  y,  fixed  (thinking  of  y^  as  the  majority  species). 


67 


Equation  (3.16)  then  transforms  to 


H"  + AnH'  - H = 0 (3.20) 


where  we  have  defined  H(n;A)  = H(5;a)  and  g(ri;A)  = g(5;a).  Note  that 


g(ri;  A) 


i H'(h;A)  ^ 

H(n;A) 


We  are  now  in  a somewhat  better  position  to  investigate  the 
asymptotic  limits  in  A of  the  behavior  of  g.  Referring  to  Eq.  (3.20) 
it  follows  at  once  that  H(q;A)  = cosh  q (from  which  we  obtain 
g(n;A)  = tanhq/n)  in  the  limit  A 0,  in  agreement  with  Eq.  (3.4).  In 

the  limit  A -*•  «>  the  gate  limit  at  n ~ 0 is  retained  by  the  factor  in 
the  middle  term  in  Eq.  (3.20).  We  consider  the  asymptotic  behavior  at 
the  other  end,  namely  as  A ->■  ».  Making  the  Ansatz,  p<2,  it 

follows  that  Eq.  (3.20)  takes  on  the  asymptotic  form 


AnH’  - H = 0 (3.21) 

as  n + from  which  we  obtain  p = 1/A,  implying  g(n;A)  » 1/Aq^. 

Evaluating  g(n;A)  by  means  of  the  power  series  expansion  for  a 
representative  sequence  of  values  of  A,  the  results  are  graphed  in  Fig. 
3.3.  The  asymptotic  limit  as  A 0 is  seen  to  carry  down  to  at  least 

two  orders  of  magnitude  in  g for  A < .01.  The  asymptotic  limit  as 
A -»■  « is  seen  to  hold  in  all  applicable  cases  ( A > 0.5).  In  the 
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o 


Figure  3.3.  g(n;A)  vs.  n for  A = 0 - 0.01,  0.1,  1.0,  10,  100,  1000  and  10 
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Intermediate  region  (0  < 4 < 0.5)  the  limits  of  precision  of  the 

calculation  suggest  that  this  limit  may  also  hold  here  as  well  but  do 
not  go  far  enough  to  confirm  it. 


We  may  define  a dimensionless  reaction  time,  most  conveniently 

as  the  point  of  intersection  of  the  as)rmptotes  in  g.  For 

A < 1,  Pc  ~ 1.  For  A > 1,  =>  1//A.  Due  to  the  effect  of  smooth 

transition,  the  prescription  in  terms  of  intersection  of  aymptotes  is 
defaulted  in  the  transition  interval  in  p.  The  retreat  of  p^  to  zero 
as  A increases  is  in  agreement  with  physical  intuition-namely , that 
for  a fixed  concentration  of  scavengee,  as  the  scavenger  concentration 
increases,  the  quickness  with  which  the  scavengee  is  removed  increases. 

Finally,  we  may  state  these  results  in  conventional  units: 


“ l/>^ky2  , 0 < A < 1 (y2  < y^  < Zy^) 

(3.22) 


Diffusion 

In  this  section  basic  diffusion  results  are  reviewed  in  order  to 
identify  the  role  displayed  by  diffusion  in  diffusion  kinetics,  and  also 
to  provide  an  independent  check  on  the  accuracy  of  the  diffusion 
component  of  the  computer  program,  DIFKIN.  More  about  the  program 
Itself  is  to  be  found  in  Appendix  A. 

Only  the  case  of  spherical  geometry  is  considered.  For  higher 
incident  energies  where  more  elongated  contours  obtain,  the  spherical 
approximation  becomes  less  adequate  and  a full 


two  dimensional 
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calculation  Is  desirable.  Efforts  to  carry  out  such  a program,  which 
are  unfinished,  are  described  In  Appendix  B.  At  lower  Incident 
energies,  the  contours  are  less  elongated,  and  the  spherical 
approximation  Is  adequate. 

For  the  sake  of  application,  an  analytic  approximation  Is  devised 
for  the  time  dependent  solution  to  the  diffusion  equation  In  spherical 
geometry.  Finally,  the  characteristic  time  for  convergence  of  the 
solution  to  steady  state  Is  discussed  In  the  context  of  providing  an 

approximate  analytic  representation  of  the  Inverse  solution  near  In  the 
approach  to  steady  state. 

The  diffusion  equation  is  given  by 


(t,t)  DV  c(r,t)  + 47rDp(r,t)  (3.23 

where  4rDp(r,t)  Is  the  source  tern.  «e  win  condense  this  Into  the 
single  dependent  variable  y(?,t).  Adequate  derivations  of  the 
diffusion  equation  In  terms  of  differential  elements  of  mass  flux  are 
given  in  various  texts,  such  as  Relf  (1965). 

Morse  and  Feshbach  (1953)  give  the  solution  to  the  aero  initial 
value  problem  by 

- I!  g„(R,r)p(?_,,t)  d?_^dt^  (3.24J 

where  g|j(8,T)  is  the  associated  Green's  function,  given  by 
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4ttD 

(2/ttDt)^ 


/4Dt 


0(t) 


(3.25) 


where  R = I I r - r 


, T = t - t . 

o 


For  the  spherical  case  we  may  set 
(3.24)  simplifies  to 


2 o’ 

F = /r  - 2rr^p^  + r^  where 

Integrating  over  y and  cb  , Eq. 

o o ^ 


t dt  00  ^ 

(r,t)  = / — ^ / dr  y(r  ,t  ) -^  fe 

0 /ttDt  0 o o r L 


(r-r  )‘ 
o 

4Dt 


4Dt 

- e“  ] (3.26) 


we  set  y(r,t)  - Ae-»"f(t)  . (6\/8,)e-6'f(t) , In  n„e  with  the  spatial 
dependence  of  radiological  yields  as  determined  In  Chapter  2.  For  this 
spatial  dependence  of  y(r,t),  the  Integral  over  r^  Is  analytic,  and 
c(r,t)  can  be  found  by  the  single  time  Integral, 


. t 

c(r,t)  = — / dt^  f(t^)[(20DT  + r)e^^^°’^  ^^^erf c 

- (26Dt  - r)e^^®°''  " ^^erfc  ~ ^)]  (3  37) 

/4Dt  ’ ^ 

The  first  term  in  the  Integrand  is  unusable  for  numerical  calculation  in 
its  present  form  for  S(0Dx  + r)  > 174,  due  to  overflow  in  the  argument 
of  the  exponent,  the  second  term  for  6(6Dt  - r)  > 174.  In  either 
case,  the  asymptotic  expansion  (Dwight,  1961)  must  be  used  to  overcome 
this  problem  - the  implied  product  of  a very  large  and  a very  small 

number  then  becomes  explicit  and  cancels  out.  These  two  terms  are  then 
given  by 
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(2DBt  + r) 


^B(6Dt  + r) 


erfc 


23Dt  + r 
/4Dt 


j-r^/4DT  s|-2B^-4  rj 

/4Dt 


(3.28a) 


and 


(2D6t-  r)  erfc  ~ " - 

/4Dt 

sgn(23DT  - r)  /4Dt/tt  e~^  ~ — ~) 

/4Dt 

+ (1  - sgn(26DT  - r))(2SDt  - r)  (3.28b) 


where 


s(x)  - z MJMi 

n-O  (Zx)^” 


Implementing  these  replacement  forms  where  needed,  the  Integral  can  be 
calculated  as  a single  integral  numerically  everywhere  (in  (r,t)). 

For  the  initial  pulse,  f(t)  = 6(t),  and  the  solution  is  analytic. 
For  the  continuous  pulse  turned  on  at  t = 0,  f(t)  = 0(t)  and  the 
solution  reduces  to  a single  integral  in  time.  In  these  calculations 
the  variable  mesh  multiple  integration  program  of  Coldweli  (private 
communication,  1979)  has  been  used.  This  program  is  relatively  cheap, 
quite  accurate,  provides  error  estimates,  and  is  convenient  to  use. 

Calculating  c(r,t)  for  the  continuous  pulse  via  DIFKIN,  agreement 
is  obtained  at  all  points  r and  t to  within  a percent.  The  solution 
c(r,t)  for  the  continous  pulse  is  graphed  in  Fig.  3.4  as  a time  sequence 


-13 


73 


o 

CO 


o 

CN 


o 

o 


o 

a\ 


o 

CO 


o ^ 
o 
u 
o 

VO 


o 

m 


o 

-<r 


o 

CO 


o 


(00 /a-fotii)  uoTtiBarju-oouoo 


Figure  3.4,  Time  sequence  of  radial  profiles  for  pure  diffusion  in  spherical  geometry. 
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of  radial  profiles.  This  sequence  graphically  displays  the  convergence 
upon  the  steady  state  profile,  Cg^(r),  which  is  found  by  repeated 
quadrature : 


A ,2(l-e'®'')  -6r, 


(3.29) 


It  is  useful  to  define  a reduced  distance,  r'  = gr  and  a reduced 
2 

time,  t'  = Dg  t.  We  may  also  introduce  a reduced  concentration,  c', 
by  the  definition,  c'  = c/(A/Dg)^,  so  that  steady  state  peak  value  is 

1. 

Regarding  the  steady  state  solution,  we  note  that  the  distance  at 
which  its  asymptotic  form  c^^(r’)  « 2/r'  equals  the  peak  value  of  the 
exact  solution,  is  twice  the  e-folding  interval.  The  inflection  point 
occurs  further  out,  at  3.38  times  the  e-folding  interval,  followed  by 
the  half  width  at  3.72.  These  distances  delineate  the  diffuse  boundary 
of  what  we  may  call  the  steady  state  core  to  reflect  the  fact  that 
steady  state  is  converged  on  fully  two  orders  of  magnitude  earlier 
within  it  than  in  the  surrounding  limb.  This  is  clearly  indicated  in 
Fig.  3.5. 

The  relative  suddenness  of  arrival  of  steady  state  within  the 
central  core  reflects  the  simple  fact  that  this  is  the  location  of  the 
source.  The  core  is  essentially  coterminous  with  the  range  of  the 
source  at  r'  = 3.5  the  source  value  is  down  to  3%  of  peak. 

It  is  possible  to  analytically  approximate  c'(r',t')  fairly  closely 
for  all  r'  and  t'  by  observing  the  early  asymptotic  limit. 


-r 


t'e  for  (r 


,2 


r'^)/4t’ 

o 


» 1 


c’(r',t') 


(3.30) 
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Figure  3.5.  Time  course  at  selected  radial  distances  for  pure  diffusion  in  spherical 
geometry.  Isocontours  of  constant  p are  overlaid. 
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where  r'  is  of  the  order  of  unity  and,  of  course  Eq.  (4.29).  Splining 
these  limits  together  by  means  of  a modification  of  the  Green's 
function,  we  arrive  at 


c'(r',t')  = c^^(r')c^(t')F  + (1  - F)[(l-G)t"  + gc^(t')]e 

where  F = F(r',t")  and  G = G(r',t") 

and  F(r',t”)  = G(r',t")  = 

where  t"  = t'/(l  + a,(4t'/(r'^  + r'^))  ^1 

i O '* 


(3.31) 


where  a - 0.56,  a = 0.65  and  r'  = 1.40  and  where  c'(t)  = c'(0,t')  is 
1 z o o 

easily  represented  by 

c^(t')  = t'/(l  + 1.40  + t').  (3.32) 

The  accuracy  of  this  fit  to  c'(r',t')  is  Illustrated  in  Figs.  3.16a  and 
b where  it  has  been  applied  to  the  time  dependent  concentration  of  H and 
OH  radical  in  the  late  chemistry  case  where  the  effect  of  chemistry  on 
initial  yield  species  is  nil.  Its  usefulness  is  further  illustrated  in 
Fig.  3.16c  (for  discussion  see  associated  text). 

In  order  to  get  a better  idea  of  convergence  time,  c(r,t)  is 
regraphed  in  Fig.  3.5  as  a family  of  time  courses  at  selected  radial 
distances  on  which  isocontours  of  convergence  to  fraction  p of  steady 
state  value  have  been  overlaid.  Bearing  the  asymptotic  form  of  the 
solution  (Eq.  3.28)  in  mind,  we  approximately  analytically  characterize 
these  isocontours  by 
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t'(r',p)  » a^(p)r'^/4  + a^Cp)  (3.33) 

where  a^^(p)  « (1  - p)^  and  a2(p)  « 1.2(1  - ( 1-p) ^ * ^ )/( 1-p) 

This  approximate  characterization  of  the  inverse  solution, 

t'(r',p),  in  the  convergence  limit  progressively  loses  accuracy  for  p < 

.25,  and  is  not  intended  to  be  correct  for  p > .95,  although  it  may  be 

good  somewhat  beyond  that.  The  asymptote  intersection  point  for  r'  = 0 

occurs  at  t’  = 1.  From  either  Eq.  (3.32)  or  (3.33)  it  follows  that 

p'  =*  0.3  at  this  point,  i.e.  at  t'  = 1.  This  may  serve  as  a suitable 

definition  of  diffusion  convergence  time,  which  we  will  denote  by  x . 

d 


Diffusion  Kinetics 

In  this  section  we  consider  the  diffusion  kinetics  of  very  simple 
systems  of  one  and  two  body  decay  of  one  species  (in  an  inert  medium) 
and  of  two  body  decay  of  two  species.  Possible  reaction  contexts  in 
which  these  kinetic  modes  are  applicable  are 


->  t 


(3.33.1) 


or 


+ M 


S2  + 


M + 


(3.33. l.M) 


and 


Si  + S2 


->  + 


k 


(3.33.2) 
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or 


+ $2  + M — ^ > S3  + M + ...  (3.33.2.M) 

respectively.  In  the  second  and  third  body  cases  (catalysis  by  M)  the 
effective  one  or  two  body  reaction  rate  coefficient,  respectively,  is 
given  by  k = C^^k'. 

The  corresponding  diffusion  kinetic  equations  are  given 
respectively  by 


8c^/at  - D^V^c^  - kcj  + y^  (3.34.1.1) 

ac^/9t  = D^7^c^  + kc^  + y^,  1 > 2 (3.34.1.2) 


and 


ac^/at  = D^V^c^  - kc^c^  + y^  (3.34.2.1) 

ac2/at  = D^V^C^  - kc^c^  + y^  (3.34.2.2) 

ac^/at  = D^v^c^  + kc^c^  + y^,  i > 3 (3.34.2.3) 

where,  in  Equations  (3.34.1.2)  and  (3.34.2.3)  it  has  been  assumed  that 

the  product  species  are  inert  and  do  not  react  further. 

Only  the  case  of  continuous  pulse  source  terms  turned  on  at  t=0  is 

considered.  In  both  single  body  and  two  body  cases  we  take 
““3r 

y^(t)  — • We  may  start  with  the  single  body  case  first. 
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Single  Body  Decay 

For  comparison  we  recall  the  time  dependence  for  the  pure  kinetic 
case  given  by  Eq.  (3.3).  The  appropriate  point  at  which  to  make  this 
comparison  is  r=0.  We  may  table  this  result  momentarily,  and  next  note 
the  spatial  dependence  of  the  steady  state  solution,  for  which  Prof. 
Green  has  provided  the  derivation. 

Setting  = C = G/r,  B = A/D  and  / = k/D,  Eq.  (3.34.1.1)  in  the 
steady  state  transforms  into 


G"  - k^G  + rBe  = 0 


(3.35) 


for  which  the  solution  is  given  by 


G 


2BB 


(e^-<^)^ 


[ 


-Kr 

e 


Bre 


-Sr 


wherefore 


(3.36) 


C 


2B3 


2 2 2 
(6  -K  ) r 


[ 


-Kr 

e 


B e 


-6r 


2 

6 -K 


(3.37) 


Note  the  reduction  to  the  case  of  pure  diffusion  for  k = 0.  The 
critical  case  3 = k is  found  by  taking  the  limit  S-k  ->■  0,  and  is  given 
by 


B e 

48' 


-3r 


(1  + 3r) 


(3.38) 


We  consider  the  following  three  cases:  (1)  3 = .25  cm  , 

.25  cm  , (2)  3 = .25  cm  \ K = .50  cm  \ and  (3)  3 = .50  cm~\ 


K 


K 
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” X o 

K = • 25  cm  . The  value  of  D used  is  1827  cm  /sec,  the  binary 

diffusion  coefficient  of  OH  in  H2O  at  T = 400  K for  medium  density 

6x10^^  molec/cc.  The  values  of  k implied  then  by  k and  D are  114.19 

sec  and  456.75  sec  , equal  to  Dg  and  4D3  , respectively,  thus 

equal  to  1/x^  and  respectively,  (see  preceeding  section).  In 

Fig.  3.6  the  time  dependence  at  the  point  r=0  as  found  via  DIFKIN  is 

plotted  for  these  three  cases  along  with  the  pure  kinetic  case  for 

comparison.  The  difference  between  the  results  for  cases  (2)  and  (3) 

which  should  and  do  approach  the  same  steady  state  limit  is 

interesting.  Agreement  to  this  point  is  to  within  a few  percent. 

In  Fig.  3.7  the  steady  state  profiles  are  plotted  along  with 
numerical  profiles  at  t = 2 sec,  2 sec  and  0.5  sec,  the  last  time 
calculated  for  cases  (1),  (2)  and  (3),  respectively.  The  agreement  is 
almost  exact.  It  is  interesting  to  note  the  exponential  fall  off  of 
these  steady  state  profiles.  This  reflects  the  independence  of  the 
logarithmic  rate  of  single  body  decay  on  concentration  while  the 
logarithmic  rate  of  diffusion  (the  so-called  "diffusion  velocity"^)  is 
independent  of  concentration  (and  can  therefore  balance  the  logarithmic 
decay  rate)  only  for  exponential  fall  off. 


^Diffusion 


"velocity" 


is  defined  by 


V. 

1 


9C. 

X 


9r 


where  D.  is  the  binary  diffusion  coefficient  of  species  i through  the 
medium  (Kanury,  1975).  The  significance  is  that  D is  defined  as 
negative  the  ratio  of  particle  flux,  F,  to  concentration  gradient. 
Therefore  F = CV. 


mole  fraction 
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Figure  3.6.  Time  course  for  single  body  decay  at  r = 0. 

Three  cases  and  pure  diffusion  shown  for 
comparison. 


mole  fraction 
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r (cm) 

Figure  3.7.  Steady  state  profiles  for  single  body  decay. 
Three  cases  as  in  Fig.  3.6. 
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In  Figs.  3.8-10  time  sequences  of  the  radial  profiles  have  been 
plotted  for  all  three  cases.  Again,  noting  the  Independence  of 
logarithmic  rate  of  single  body  decay  on  concentration,  we  observe  In 
the  two  cases  where  the  chemistry  Is  as  fast  as  the  diffusion  at  the 
center  (k  > 8)  that  the  convergence  Is  essentially  simultaneous  at  all 
spatial  points.  In  the  case  In  which  the  chemistry  Is  slower  than  the 
diffusion  at  the  center,  this  Is  seen  not  to  be  the  case.  This  Is 
because  the  early  time  e fall  off  Is  relaxing  to  the  later  time 
e /r  fall  off.  The  trend  towards  the  pure  diffusion  limit  Is  obvious. 
In  the  fast  chemistry  case,  the  early  time  e”^^  fall  off  remains 
seated.  With  a little  effort  a suitable  analytic  spline  could  be  found 
to  patch  early  and  later  time  dependences  together,  and  thereby 
approximate  the  time  dependent  solution,  c(r,t). 

Two  Body  Decay  of  One  Species 

We  next  turn  our  attention  to  the  two  body  case.  We  first  consider 
the  case  of  two  body  decay  of  one  species  with  Itself.  To  reflect 
counting  (stoichiometry)  the  destruction  term  Is  multiplied  by  two: 

-^  = DV  C - 2kC^  + y (3.39) 

We  again  use  D = 1827  cm  /sec  and  3 = .25  cm  and  we  use  the 
following  four  values  of  k:  1.08x10^^  cc/mole/sec,  6.31x10^^ 

cc/raole/sec,  1.08x10^^  cc/mole/sec  and  1.08x10^^  cc/mole/sec.  The  value 
of  A Is  1.03x10  mole/cc/sec.  The  reaction  times  for  pure  kinetics  at 
r=0  Implied  by  these  numbers,  using  Eq.  (4.4)  are  .021  sec,  .0088  sec. 
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. Time  sequence  of  radial  profiles  for  single  body 
decay,  6 = .25  cm”^,  k = .25  cm~^. 


Figure  3.8 
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Figure  3.9.  Time  sequence  of  radial  profiles  for  single  body 
decay,  3 = . 25  cm~^,  k = . 50  cm“  . 
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. Time  sequence  of  radial  profiles  for  single  body 
decay,  6 = .50cm“l,  k = . 25  cm“l. 


Figure  3 . 10 
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.0021  sec  and  .00021  sec,  respectively.  The  diffusion  convergence 
2 

time  Tjj  = 1/D0  = .0088  sec.  These  four  cases  then  give  one  slow 

chemistry,  one  with  chemical  reaction  time  equal  to  diffusion 
convergence  time,  and  two  with  fast  chemistry. 

Making  the  transformation  C = G/r,  and  considering  only  the  steady 
state  limit,  Eq.  (3.39)  becomes 


G" 


2 

2k  9_ 

D r 


Bre 


-Br 


= 0 


(3.40) 


We  may  define  an  Inverse  length  k associated  with  k/D  but  in  a 

different  way  than  in  the  single  body  case.  We  may  ask  that  be 
4 

replaced  by  k reflecting  the  number  of  bodies.  We  may  allow  it  to 
become  a coefficient  of  the  second  term  by  admitting  a compensating 
coefficient.  Finally,  we  may  look  somewhere  else  for  analogies — namely, 
we  may  compare  the  pure  diffusion  steady  state  peak  (A/D6^  = B/B^)  with 
the  pure  kinetic  steady  state  limit  ( A/2k).  This  suggests  that  we 
define  k = ‘^/2kA/D^  = ^/2kB/D  so  that  B/k^  = /A/ 2k.  Equation  (3.40) 
then  becomes 


G" 


■4^  + - 0 


(3.41) 


This  is  all  very  well  and  good,  but  what  has  it  gained  us?  We  look 
at  Fig.  3.11  where  the  time  dependence  of  the  solution  at  r=0  is  plotted 
for  the  four  cases  listed  above  and  for  pure  diffusion  for  comparison. 
For  the  four  values  of  k listed  above,  k = .16,  .25,  .51  and  1.6.  Here 
now  is  the  rest  of  the  explanation  for  the  choice  of  the  values  of  k: 
we  find  that  by  defining  k as  we  did,  so  that  the  forms  for  steady 


n)Qle  fractioa 
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t(sec) 


Figure  3.11.  Time  course  for  two  body  single  species 
decay  at  r = 0.  Three  cases  sho\<m. 
Corresponding  pure  kinetic  cases  and  pure 
diffusion  shown  for  comparison. 
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state  pure  diffusion  peak  value  and  steady  state  pure  kinetic  steady 

state  peak  value  are  analogous,  a 'characteristic  diffusion  time', 

2 

1/D3  is  indicated  by  analogy:  8 = < =>  Noting  the  values 

of  the  steady  state  limits  at  r=0  for  the  assorted  diffusion  kinetic 

cases,  we  find  that,  to  an  average  accuracy  of  < 10%,  this  limit  is 
2 2 

given  by  C = B/(3  + k ).  We  may  refer  to  the  case  k « 8 (t  » t ) 

c d 

as  weak,  slow  or  late  chemistry,  and  to  the  case  ic  » 8 (t  « t ) as 

c d 

strong,  fast  or  early  chemistry. 

The  time  evolutions  of  profiles  for  the  first,  third  and  fourth 
cases  listed  above  are  plotted  in  Figs.  3.12-14.  The  steady  state 
profiles  of  all  four  cases  along  with  the  steady  state  profile  for  pure 
sre  plotted  in  Fig.  3.15.  A striking  and  very  interesting 
array  of  patterns  meets  us  here.  Looking  at  the  time  sequences  of 
profiles,  we  observe  that  as  the  chemistry  becomes  stronger,  the 
evolution  of  the  steady  state  profile  resembles  the  stripping  away  from 
an  Irreducible  core  of  more  and  more  covering. 

What  is  happening  is  this.  In  the  core,  as  the  chemistry  becomes 
faster,  it  seats  a steady  state  profile  long  before  diffusion  can  have 
an  effect,  resulting  in  a purely  local  chemical  effect.  Thus  in  the 
core  the  effect  of  the  diffusion  term  drops  off,  leaving  a local  balance 
between  destruction  term  and  source  term. 

In  the  limb  the  effect  of  the  short  range  source  term  dies  away. 
Without  chemistry  the  pure  diffusion  1/r  limit  would  take  place, 
reflecting  only  diffusion.  With  chemistry,  a balance  is  set  up  between 
diffusion  and  chemical  destruction: 


2 

2k  G2 
D r 


G 


= 0 or  G" 


0 


(3.42) 
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Figure  3.12.  Time  sequence  of  radial  profiles  for  two  body  single  species  decay, 
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Figure  3. 13.  Time  sequence  of  radial  profiles  for  two  body  single  species  decay. 
6 = .25  cm~^.  K .51  cm"'. 
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Figure  3.14.  Time  sequence  of  radial  profiles  for  two  body  single  species  decay. 
3 = .25  cm~^ . K = 1.6  cm""!. 


93 


r Ccm) 


Figure  3.15. 


Steady  state  profiles  for  two  body  single  ^ 
species  decay:  3 = . 25  cm~^.  k = . 16  cm  , 

.25  cm~^,  .51  cm“l  and  1.6  cm~l.  Overlaid 
by  fit  obtained  in  Appendix  D. 


94 


This  asymptotically  limiting  form  of  the  equation  is  solved  by  the 
power  law  form,  G = 2B/k\  = D/kr  (wherefore  C = 2B/<\^  = D/kr^), 
which  is  independent  of  A,  the  scaling  of  the  source.  This  invariance 
is  a result  of  the  nonlinear  dependence  on  concentration  of  the 
diffusion-kinetic  balance.  Looking  to  Fig.  3.15  we  note  that  the 
solution  as  calculated  numerically  is  indeed  approaching  the  asymptote, 
except  for  the  weakest  chemistry  case,  in  which  a 50%  discrepancy 
appears.  This  has  not  yet  been  explained  and  may  be  a numerical 
artifact. 

Proceeding  from  these  asymptotic  trends  in  k/B  and  r,  it  is 
possible  to  construct  an  analytic  form  which  splines  them  together  over 
the  entire  range  of  ic/S*  This  has  been  done,  and  the  derivation  and 
form  are  presented  in  Appendix  D.  Resultant  values  are  overlaid  on  the 
true  values  in  Fig.  3.15.  Agreement  is  very  good.  With  a little  effort 
it  should  also  prove  possible  to  find  a suitable  analytic  spline  which 
patches  together  early  and  later  time  asymptotic  limits. 
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Two  Species  Two  Body  Decay — Specific  Example 
In  the  case  of  two  body  decay  coupling  two  distinct  species,  two 
new  degrees  of  freedom  are  introduced:  (1)  y^^  t y^,  and  (2)  * D^. 

We  consider  the  reaction 

H + OH  + M -*■  H^O  + M (3.43) 

with  T = 400  K,  n = 1.84x10^^  #/cc  (p  = 1 atm).  The  source  terms,  y^^ 

and  yQjj,  are  taken  to  be  the  post-recombination  spatial  yields  obtaining 

for  a 1 keV  electron  beam.  These  stand  in  the  ratio,  y^/y2  = The 

diffusion  coefficients,  found  from  Table  6.1,  stand  in  the  ratio, 

Dj/Dz  = 3.5.  Curiously  enough,  as  we  shall  see  below,  the  more  powerful 

effect  is  from  the  differential  in  source  term. 

Three  injection  rates  are  considered:  I = = 3.67x10^*^  ///sec, 

I = 100  Ip  and  I = 10^  x Ip  Reaction  times  are  1.3xl0~^  sec,  1.3xl0~^ 

sec,  and  1.3x10  ^ sec,  respectively.  Diffusion  convergence  time  is 
— ft 

1.3x10  sec.  Thus,  x = lOx.,  x = x,,  and  x = x,/10,  or 

c a c d c d 

2 2 2 2 2 2 

equivalently,  g = 10<  , 6 = < and  g = K /lO,  respectively.  We 

thus  have  one  representative  case  each  of  slow,  borderline  and  fast 
chemistry.  The  problem  is  abstracted  from  the  three  cases  at  one 
atmosphere  pressure,  T = 400  K,  and  for  = 1 keV  considered  in  Chapter 
VI.  Although  we  did  not  consider  the  evolution  of  the  concentration  of 
product  species  in  the  previous  section,  we  do  so  here,  labelling 
product  H2O  separately  from  the  medium. 

Time  course  of  the  system  at  r = 0 is  plotted  for  these  three  cases 
in  Figs.  3.16-18.  Time  sequences  of  radial  profiles  are  plotted  in 
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Figure  3.16.  Time  course  at  r = 0 for  H + OH  + M ^ H 0 + M. 

Slow  chemistry  case.  Compare  Figure  6.^. 
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Figure  3.17.  Time  course  at  r = 0 for  H + OH  + M ->  H2O  + M. 

Borderline  chemistry  case.  Compare  Figure  6.8. 
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Figs.  3.19a-c,  3.20a-c,  and  3.21a-c.  Let  us  consider  the  first  case,  the 
slow  chemistry  case  (see  the  first  time  plot  and  the  first  group  of 
radial  profile  time  sequences). 


The  Slow  Chemistry  Case 

In  the  slow  chemistry  case  diffusion  shaped  profiles  have  converged 
before  reaction  time  at  center  is  reached.  From  this  point  on,  It 
becomes  relevant  to  compare  with  the  reaction  time  for  initial 
concentrations,  thinking  in  terms  of  reinitialized  time.  The  reaction 
time  follows  from  Eqs.  3.6  and  3.9  and  is  given  by  Eq.  (4.11),  and, 
substituting  values  of  concentrations  and  reaction  rate  coefficient,  we 
obtain  = 2.4  x 10  sec.  Why  then  is  product  H2O  concentration  at 
r = 0 is  still  not  completely  converged  at  t = 10“^  sec?  There  are  two 
reasons:  (1)  convergence  is  more  gradual  in  the  initial  pulse  case 

(compare  Figs.  3.1  and  3.2),  and  (2)  the  source  term  of  product  H2O 
continues  to  spread  outward  after  it  has  converged  in  the  center.  H2O 
produced  in  the  limb  diffuses  towards  the  center  as  well  as  outward.  As 
a result  the  product  H2O  concentration  at  the  center  continues  to 
increase,  prolonging  the  convergence. 

Since  the  time  courses  of  the  reactant  concentrations  are  those  for 
pure  diffusion,  for  which  we  have  an  accurate  analytic  representation, 
we  may  calculate  the  time  course  of  the  product  concentration  by  means 
of  the  Green's  function  integral.  Results  are  overlaid  in  Fig.  3.19c. 

(To  illustrate  the  quality  of  the  fit  to  the  diffusion  time  courses,  the 
pure  diffusion  fit  is  overlaid  in  Figs.  3.19a  and  b.  The  agreement  of 
the  calculated  product  concentration  predicted  on  the  basis  of  the 


Figure  3.19.  Time  sequence  of  radial  profiles  for  (a)  H,  (b)  OH  and 
(c)  H2O  in  the  slow  chemistry  case.  Overlaid  by  pure 
diffusion  fit  in  (a)  and  (b)  and  by  Green's  function 
integral  calculated  value  in  (c)  based  on  pure 
diffusion  fits  in  (a)  and  (b). 
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Figure  3.19  (c) 
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Figure  3.19  (c ' ) 
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Figure  3.20.  Time  sequence  of  radial  profiles  for  (a)  II,  (b)  OH. 
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Figure  3.20.  Time  sequence  of  radial  profiles 
for  (c)  H„)  in  the  borderline 
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fitted  reactant  concentrations  is  remarkably  good.  With  it  we  can 

predict  the  time  course  out  to  very  distant  times  out  to  which  it  would 

be  impractical  to  calculate  using  the  marching  step  code.  Results  at 

these  later  times  are  also  plotted  in  Fig.  3.19c. 

Looking  at  the  overlay  in  Fig.  3.19c,  or  in  Fig.  3.19c’  where  it  is 

taken  out  to  much  longer  times,  it  becomes  apparent,  assuming  this 

numerical  integral  is  correct,  that  the  growth  becomes  logarithmic. 

This  result  is  easily  confirmed  analytically.  The  source  term  is 

proportional  to  1/r^  in  the  limbs,  cancelling  the  volume  element  scale 

factor.  The  effective  upper  limit  for  r is  /4Dx.  Partial 

3/2 

cancellation  of  the  (4Dt)  in  the  denominator  of  the  x-integral 
then  yields  the  logarithmic  result.  The  numerical  result  is,  therefore, 
correct. 


Borderline  and  Fast  Chemistry  Cases 

We  next  consider  the  two  other  cases  of  borderline  chemistry  and 
fast  chemistry.  In  contrast  with  the  sequence  of  cases  in  the  single 
species  case  (previous  section)  where  the  reaction  rate  coefficient  was 
varied  while  the  source  term,  y,  remained  fixed,  here— in  a physically 
realistic  context  the  rate  coefficient  remains  fixed  (for  fixed  density 
the  effective  two  body  rate  coefficient  implied  by  the  three  body 
coefficient  Is  fixed),  while  the  source  terra  is  varied.  Comparison  is 
therefore  aided  by  scaling  concentrations  to  source  terra.  In  this 
context  it  is  relevant  to  consider  Fig.  3.11  where  diffusion  kinetic 
time  courses  at  r=0  are  compared  with  the  corresponding  pure  kinetic 
courses  and  the  corresponding  pure  diffusion  time  course  for  the  single 


species  case. 
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Although  it  would  be  informative  to  similarly  graph  corresponding 
time  courses  together  on  one  graph  here,  we  can  derive  the  essence  of 
this  information  from  Table  3. 1 in  which  steady  state  peak  values  are 
compared. 

In  the  borderline  case,  the  peak  levels  of  both  reactant  species 
are  diminished  by  the  same  relative  amounts.  Both  profiles  exhibit  the 
stronger  convergence  to  steady  state  near  the  peak  characteristic  of 
increasing  chemical  effect.  The  fall  offs  of  both  profiles  have  become 
steeper  by  a small  amount,  somewhat  more  for  OH  than  for  H,  even 
relative  to  its  extent. 

In  the  fast  chemistry  case,  the  OH  peak  level  is  decreased  by  twice 
as  much  as  is  the  H peak,  bringing  them  closer  to  their  source  term 
ratio,  1.1:1.  This  increased  localization  reflects  the  greater  speed  of 
the  chemistry  outpacing  the  slower  diffusion  (nonlocalizing  effect). 

The  H fali  off  becomes  even  broader,  while  the  OH  fall  off  becomes 
extremely  steep,  in  fact,  essentially  exponential  fall  off.  The  greater 
diffusivity  of  H is  not  the  source  of  this  effect.  The  differential  in 
the  source  terms  is.  Recall  the  results  for  two  species  two  body  decay 
for  unequal  source  terms  in  the  pure  kinetic  case.  In  the  steady  state 
limit  the  minority  species  decays  to  zero  while  the  majority  species 
then  grows  at  the  rate  of  the  difference  of  their  source  terms.  In  the 
diffusion  kinetic  case,  diffusion  braking  sets  a limit  on  the  scavenging 
strength  of  the  majority  species,  and,  as  a result,  a nonzero  fraction 
of  the  minority  species  escaped  destruction,  while  the  majority  species 
enjoys  a pure  diffusion  limit  in  the  limb. 
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TABLE  3.1 

Comparison  of  Steady  State  Peak  Concentrations 

(first  entry  is  mole  fraction;  second  entry  is  mole  fraction 
normalized  to  injection  rate  of  slow  chemistry  case;  third 
entry  is  this  quantity  normalized  to  slow  chemistry  mole  fraction) 


slow  chemistry 


borderline  chem. 


H 

OH 

2.0(-6) 

2.0(-6) 

1 

5.6(-6) 

5.6(-6) 

1 

1.2(-4) 

1.2(-6) 

.6 

3.6(-4) 

3.6(-6) 

.6 

4.0(-3) 

4.0(-7) 

6.0(-3) 

6.0(-7) 

.2  .1 


fast  chemistry 
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The  gradual  1/r  fall  off  of  the  majority  species  has  about  the  same 
effect  as  a constant  level  and  the  steady  state  limit  profile  of  the 
minority  species  Is  essentially  exponential,  mimicking  single-body 
decay. 


Fast  Chemistry  with  Reactant  Yield  Rates  Interchanged 

In  order  to  demonstrate  conclusively  that  the  effect  seen  here  Is 
due  to  the  10%  difference  In  the  yield  rates,  rather  than  to  the  large 
difference  In  the  diffusion  constants,  the  same  system  was  recalculated 
with  the  yield  rates  Interchanged.  The  time  plot  Is  given  In  Fig.  3.22, 
and  the  profiles  are  given  In  Fig.  3.23.  It  can  be  seen  at  once  that 
the  differential  in  the  yield  rates  Is  the  source  of  the  effect.  Here 
the  OH  enjoys  the  pure  diffusion  limit  In  the  limb,  while  the  H is  held 
down  to  an  essentially  exponential  fall  off.  The  strange  behavior  in 
the  later  time  H profiles  In  the  limb  may  be  a numerical  artifact  or 
partly  so.  This  question  could  be  resolved  by  increasing  the  number  of 
radial  grid  points  and  extending  the  edge  point  outward,  and  replacing 
the  Simpson  rule  Interpolation  by  cubic  Interpolation  (see  Appendix  B). 

We  note  that  the  H2O  profiles  are  only  minimally  perturbed  by  the 
Interchange  In  yield  rates.  This  makes  sense  since  the  H2O  production 
rate  Is  proportional  to  the  product  of  the  H and  OH  concentration 
levels.  Thus,  to  start  with,  at  early  times  before  diffusion  has  an 
effect,  nothing  is  changed  for  H2O  production.  At  later  times  it  can  be 
seen  from  the  profile  plots  that  the  bulk  of  newly  created  H and  OH  go 
into  H2O  production;  therefore,  there  Is  little  change  here  either. 


concentra,t ion.  (njole/cc) 
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Figure  3.22.  Time  course  at  r = 0 for  H + OH  + M H„0  + M. 

Fast  chemistry  case.  H and  OH  yield  rates 
reversed.  Compare  with  Figure  3.15. 


concentration  Cmole/cc) 
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•01  .02  .03 

r Ccm) 


Figure  3.23 


Time  sequence  of  radial  profiles 
for  (a)  H. 
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0 .01  . .02  .03 
r Ccm) 

Figure  3.23.  Time  sequence  of  radial  profiles 
for  (b)  OH. 
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Figure  3.23.  Time  sequence  of  radial  profiles 
for  (c)  H2O  in  the  reversed  fast 
chemistry  case.  Compare  with  Figure 
3.18. 
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We  note,  finally,  the  still  lower  level  of  H radical  at  r=0  In  the 
steady  state  limit  as  a result  of  the  Interchange  of  yield  rates. 

Whereas  It  was  lower  before  because  of  Its  greater  dlffuslvlty  sending 
It  out  to  the  limb.  It  Is  still  lower  here  as  a result  of  the  loss  of 
advantage  In  source  term.  The  hump  In  the  time  course  at  r=0  and  the 
broad  swing  out  and  back  In  the  profile  reflect  further  this  association 
of  high  dlffuslvlty  with  disadvantage  In  yield  rate.  In  the 
Intermediate  time  frame  of  this  swing,  the  more  diffusive  minority 
species  (H)  has  a greater  survival  rate  In  the  limb  where  the  more 
slowly  moving  majority  species  (OH)  has  not  yet  reached.  This  phase 
comes  to  an  end,  however.  In  a delayed  outward  sweep  of  OH  chasing  after 
this  escaped  prey.  The  pinched  off  bubble  of  H moving  outward  ahead  of 
the  OH  may  therefore  be  at  least  partly  real. 


Where  We  Go  From  Here 

This,  finally.  Is  a good  point  at  which  to  pause  for  a moment  and 
take  stock  before  proceeding.  We  have  laid  out  the  mathematical 
formulae  of  simple  kinetics  of  one  and  two  body  decay  for  both  initial 
pulse  and  continuous  steady  pulse  source  terms,  for  both  one  species  and 
two  species  reactive  systems.  We  have  calculated  the  complete  time 
dependence  of  pure  diffusion  for  spherical  geometry,  and  characterized 
It  In  sufficient  detail  to  aid  us  In  our  diffusion  kinetic  work. 

Since  the  overall  features  of  both  diffusion  and  kinetics  are 
Intuitively  straight  forward,  the  Important  benefit  of  this  work  was 
quantitative  identification  of  reaction  times  and  diffusion  convergence 


time. 
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these  basic  results  at  our  disposal,  we  then  went  on  to 
consider  diffusion  kinetics  of  one  and  two  body  decay.  More  attention 
was  paid  to  two  body  decay.  However,  effective  one  body  decay  through 
reaction  with  the  medium  (H2O)  does  occur  In  both  the  Ion  and  the 
neutral  reactive  systems. 

Two  body  decay  was  first  treated  In  the  simpler  context  of  single 
species.  Analytic  asymptotic  limits  were  identified.  The  interplay  of 
diffusion  and  kinetics  was  clearly  discerned.  The  two  species  case  was 
considered  in  the  specific  physical  context  of  the  reaction,  H + OH  + M 
H^O  + M,  for  realistic  radiological  inputs,  in  fact  for  one  of  the 
density-temperature-incident  energy  configurations  which  will  be 
considered  in  Chapter  VI,  so  that  isolation  of  effect  is  possible. 

What  remains  to  be  done  is  as  follows:  1.  The  instantaneous 

of  the  ion  reactions  must  be  put  in  place  before  we  can 
proceed  to  the  neutral  reaction  set  operating  by  itself.  2.  The 
kinetics  of  the  neutral  reaction  set  are  to  be  identified.  Since  we 
view  diffusion  kinetics — in  the  context  of  the  type  of  problem 
considered  here  as  diffusion  limitation  or  modification  of  kinetics,  we 
need  kinetics  to  modify.  Diffusion  kinetics  of  the  full  neutral 
reaction  set.  Here  we  tie  back  in  to  the  spatial  yields,  and  finally, 
into  their  energy  dependence,  and  therewith  the  energy  dependence  of  the 


diffusion  kinetics. 


CHAPTER  IV 
ION  REACTIONS 

The  ion  reactions  divide  into  three  categories:  (1)  ion— neutral 
reactions,  for  which  only  those  in  which  H2O  is  the  neutral  species  are 
of  any  consequence  due  to  concentration,  (2)  hydration  of  the  proton, 
which  as  we  will  see  below  has  negligible  stoichiometric  effect,  and  (3) 
proton-electron  recombination. 

In  his  Ph.  D.  thesis  (1983)  Rio  has  compiled  the  predominant  ion 
reactions  (and  associated  reaction  rates)  which  take  place  in  water 
vapor  subjected  to  electron  bombardment  (at  low  temperatures,  i.  e. , ca. 
200  - 450  K).  We  adopt  this  reaction  set  here.  Reaction  43  (H20‘^  + H2 
-*■  H^O  + H)  is  omitted  for  the  reason  mentioned  above.  These 
reactions  and  reaction  rates  are  given  in  Table  4.1.  The  sources  of 
these  reactions  and  reaction  rates  are:  Huntress  (1973)  and  Ferguson 

(1973)  for  the  ion-neutral  reactions;  Good_et  al.  (1970)  for  the  ion 
clustering  reactions;  and  Leu  et  al.  (1973)  for  the  recombination 
reactions. 

The  neutral  reactions  and  reaction  rates  are  taken  from  a recent 
combustion  oriented  compilation  of  Westbrook  (1981),  and  are  given  in 
Table  4.2.  These  are  the  same  as  used  by  Rio,  and  have  been  implemented 
in  previous  combustion  (and  radiolysis)  related  work  in  this  group 
(Green  et  al. , 1983;  Pamldimukkala  et  al.,  1983). 

It  is  a well  known  fact  that  the  ion  reactions  proceed  much  faster 
than  the  reactions  between  neutrals  for  typical  medium  densities  and 
injection  rates.  Taking  typical  concentrations  into  account,  a quick 
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TABLE  4.1 

Reactions  and  Reaction  Rates  for  the  Ions 


Reaction 

Reaction  Rate 
ca.  300-400  K 

Coefficient  at 
(cc-mole-sec) 

a 

H20''’  + H2O 

— > 

H30‘^  + OH 

1.14(15) 

b 

OH'*'  + H2O 

— > 

H20'’’  + OH 

3.49(14) 

c 

OH'*'  + H2O 

“> 

H30'^  + 0 

2.83(14) 

d 

H'*'  + H2O 

— > 

H20"''  + H 

4.94(15) 

e 

+ H2O  + M 

"> 

H'^(H20)2+  M 

1.23(21) 

f 

H‘''(H20)2  + H2O  + M 

— > 

H'^(H20)3+  M 

8.34(20) 

g 

H"^(H20)2  + e“ 

~> 

H + H2O  + H2O 

1.63(18) 

h 

H‘^(H20)3  + e“ 

— > 

H + H2O  + H2O  + H 

2O  2.77(18) 

i 

H3O  + e 

— > 

H + H2O 

7.83(17) 

j 

+ e 

— > 

H + H + OH 

«7. 83(17) 
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TABLE  4.2 

Reactions  and  Reaction  Rates  for  the  Neutrals 
(Reaction  rate  coefficients  are  in  cc-mole-sec  units.) 


Reaction 


If 

H + O2 

~> 

0 + OH 

lb 

H + ©2 

<— 

0 + OH 

2f 

H2  + 0 

— > 

H + OH 

2b 

H2  + 0 

<— 

H + OH 

3f 

H2O  + 0 

--> 

OH  + OH 

3b 

H2O  + 0 

<~ 

OH  + OH 

4f 

H2O  + H 

— > 

H2  + OH 

4b 

H2O  + H 

<— 

H2  + OH 

5f 

H2O2  + OH  - 

->  H2O  + 

HO2 

5b 

H2O2  + OH  <■ 

- H2O  + 

HO2 

6f 

H2O  + M 

— > 

H + OH  + 

M 

6b 

H2O  + M 

<— 

H + OH  + 

M 

7f 

H + O2  + 

• M - 

->  HO2  + 

M 

7b 

H + O2  + 

M <—  HO2  + 

M 

8f 

HO2  + 0 

— > 

OH  + O2 

8b 

HO2  + 0 

<~ 

OH  + O2 

Reaction  Rate  Coefficient  at 
T = 300  K 400  K 


1.09(2) 

4.73(12) 

1.24(5) 

6.29(12) 

1.76(6) 

2.16(7) 

9.95(7) 

5.30(8) 

1.45 

5.00(11) 

3.19(3) 

7.91(11) 

1.55(-1) 

3.87(9) 

7.73(2) 

3.36(10) 

4.88(11) 

3.64(-ll) 

1.04(12) 

3.43(-5) 

7.00(-60) 

1.57(18) 

9.00(-42) 

8.81(17) 

8.88(15) 

8.30(-19) 

5.84(15) 

1.90(-10) 

9.36(12) 

3.68(-28) 

1.42(13) 

7.54(-18) 
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9f  HO2  + H — > OH  + OH 
9b  HO2  + H <—  OH  + OH 

lOf  HO2  + H — > H2  + O2 
10b  HO2  + H <—  H2  + O2 

Ilf  HO2  + OH  — > H2O  + O2 
11b  HO2  + OH  <—  H2O  + O2 

12f  H2O2  + ©2  HO2  + HO2 

12b  H2O2  + ©2  — HO2  + HO2 

13f  H2O2  + M — > OH  + OH  + M 
13b  H2O2  + M < — OH  + OH  + M 

14f  H2O2  + H — > HO2  + H2 
14b  H2O2  + H <—  HO2  + H2 

15f  0 + H + M — > OH  + M 
15b  0 + H + M < — OH  + M 

16f  O2  + M — > 0 + 0 + M 
16b  O2  + M < — 0 + 0 + M 

17f  H2  + M — > H + H + M 
17b  H2  + M <~  H + H + M 


1.04(13) 

7.32C-17) 

2.30(13) 

1.47(-9) 

7.76(12) 

4.26(-29) 

1.04(13) 

1.43(-18) 

9.36(12) 

9.74(-40) 

1.42(13) 

2.76(-26) 

3.42(-18) 

1.87(12) 

2.00(-10) 

2.84(12) 

8.52(-17) 

4.51(18) 

1.65(-8) 

5.38(17) 

3.15(9) 

1.72(-2) 

1.51(10) 

43.8 

1.00(16) 

0. 

1.00(16) 

4.20(-40) 

0. 

9.47(14) 

7.40(-48) 

8.74(14) 

2.50(-56) 

3.02(15) 

7.66(-39) 

3.02(15) 
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glance  at  the  reaction  rates  in  Tables  4.1  and  4.2  makes  this  clear.  In 
order  to  properly  compare,  typical  concentrations  must  be  taken  into 
account.  In  the  interest  of  thoroughness  it  would  be  advantageous  to 
compare  over  a representative  grid  in  medium  density  and  injection 
rate.  This  is  done  in  the  last  section  of  this  chapter  with  follow-up 
in  the  last  section  of  the  following  chapter. 

Rather  than  just  look  at  isolated  reaction  rates,  it  is  much  more 
informative  to  calculate  the  evolution  of  the  full  kinetic  system  of 
ions  and  neutrals  in  order  to  identify  the  characteristic  trends  of  the 
ion  regime  and  of  the  separation  of  time  scales  of  the  two  systems.  For 
this  purpose  both  the  initial  pulse  case  and  the  continuous  pulse  case 
are  of  interest.  It  would  be  useful  to  perform  such  kinetic 
calculations  over  the  above  mentioned  representative  grid  in  medium 
density  and  injection  rate.  To  this  purpose  the  code  DIFKIN,  written  to 
calculate  diffusion  kinetics  or  pure  kinetics,  could  be  used.  However, 
a pre-existing  kinetic  code  written  by  Rio  (Ph.  D.  thesis,  1983)  based 
on  a the  Gear  method  (Hindmarsh,  1974)  has  already  proved  useful  to  this 
purpose,  as  demonstrated  by  the  initial  pulse  calculation  presented  in 
Rio  (Ph.  D.  thesis,  1983,  Fig.  3.1)  and  reproduced  here  in  Fig.  4.1. 

This  kinetic  plot  is  extremely  informative  and  we  shall  now  discuss  it. 

The  H2O  concentration  is  1 x lo“^  moles/cc.  The  Initial 
concentrations  of  the  radiological  deposition  products  are  apparent  from 
the  plot.  As  stated  by  Rio,  these  concentrations  were  chosen  to  be 
proportional  to  the  Initial  radiological  yields.  These  proportions  are 
essentially  constant  with  respect  to  as  a result  of  the  essentially 
constancy  of  the  g-values  (see  Chapter  II). 
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Figure  4.1. 


Kinetic  time  course  of  full  ion-neutral  system  for 
initial  pulse,  (Rio,  Ph.D.  thesis,  1983). 
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Inspection  of  this  plot  reveals  that  proton  transfer  from  the 

dominant  primary  ionic  species,  to  H2O,  is  well  underway  by  lO"®  sec. 
Specifically,  the  H2O  concentration  has  fallen  by  a factor  of  about  1/e 
at  this  time.  peaks  shortly  after  that,  with  H‘''(H20)2  peaking 

shortly  after  it,  and  H'^(H20)3  shortly  after  it,  spread  over  the 

interval  2 x 10  sec  < t < 2 x 10  ^ sec.  Recombination  gets  underway 

~8 

at  about  2 x 10  sec  and  is  essentially  complete  by  5 x 10~^  sec. 

The  neutral  reactions  do  not  become  critical  until  after  1 x 10~^  sec. 
Thus,  the  ion  reactions  are  well  separated  from  the  neutral  reactions  in 
time~~by  more  than  three,  almost  four  orders  of  magnitude. 

It  may  be  noted  in  passing  that  electron  recombination  with  the 
primary  yield  ions  occurs,  but  at  a rate  negligible  in  comparison  with 
the  rate  of  their  reaction  with  the  much  more  numerous  H2O  molecules. 

Inferring  from  the  comfortable  separation  of  time  scales  between 
the  ion  reactions  and  the  neutral  reactions,  we  may  conclude  that  for  a 
wide  range  of  medium  densities  and  injection  rates  that  this  separation 
would  maintain.  Where  it  may  begin  to  break  down  should  be  at  higher 
densities,  where  the  enhanced  reaction  rate  of  the  three  body  radical 
recombination  reactions  would  begin  to  push  the  onset  of  neutral 
activity  backward  in  time.  Such  a calculation  at  liquid  density,  for 
example,  should  be  made  to  span  the  interval  so  that  we  may  speak  with 
some  confidence  of  this. 


Reduced  Reaction  Set 

Because  the  ion  reactions  are  all  one~way,  and  a nondepleting 
regime  is  assumed,  ruling  out  feedback  loops  via  reactions  of  the  ions 
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with  neutrals  other  than  H2O,  the  ion  reactions  constitute  a feedforward 
system.  Reactions  b and  d (electron  transfer)  feed  into  reaction  c 
(proton  transfer)  which  parallels  reaction  a (proton  transfer).  The  end 

point  of  these  reactions  of  the  ions  with  H2O  is  the  hydronium  ion, 

H3O+. 

Noting  the  increase  in  the  recombination  reaction  rates  with 
increasing  level  of  hydration,  it  follows  that  the  effect  of  the 
multiple  hydration  (reactions  e and  f)  is  to  shorten  the  time  the  proton 
remains  imrecombined.  This  trend  in  reaction  rate,  therefore  in 
reaction  cross  section  is  apparently  due  to  a diffusion  of  positive 
charge  over  a greater  volume  via  polarization. 

The  end  step,  electron-proton  recombination  splinters  off  these 
associated  H2O  molecules,  and  finalizes  the  production  of  H.  Reaction  j 
is  energetically  unfavored,  and  so  can  be  neglected.  As  a result,  the 
multiple  hydration  reactions  have  negligible  stoichiometric  effect. 

From  the  above  analysis  we  may  predict  the  post-recombination 
concentrations  in  terms  of  the  initial  concentrations  following  an 
initial  pulse  of  the  electron  beam,  and  post-recombination  yield  rates 
in  terms  of  initial"  yield  rates  for  a continuous  beam,  assuming  the 
one  other  important  point  which  we  will  get  to  momentarily,  the 
separation  in  time  scales  of  the  ion  reactions  from  the  neutral 
reactions.  For  the  moment  this  will  be  assumed.  Let  us  first  consider 
the  case  of  an  Initial  pulse.  We  first  define  a quantity  f which 
defaults  to  1 if  the  electrons  outnumber  the  ions,  and  which  otherwise 
^^^^cts  the  electron  deficit.  Formally,  we  may  write 


f = r(l-e(r-l))  + 0(r-l) 


(4.1) 
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where 


r = (e")/[(H'^)  + (OH*^)  + (H^O*^)] 

(4.2) 

and  9 is  the  Heaviside  step  function.  The  post-recombination 
concentrations  in  terms  of  the  initial  concentrations  are  then  given  by 

[h]  = (H)  + (h"^)  + f[(OH’*')  + (h"^)  + (H20'^)] 

(4.3a) 

[H2]  = (H2) 

(4.3b) 

[oh]  = (OH)  + (2k^/(k^+k^))(0H'’')  + (h"^)  + (H20"^) 

(4.3c) 

[0]  = (0)  + (k^/(kj^+k^))(0H‘^) 

(4.3d) 

[e~]  = (e")  - f[(0H'^)  + (h'^)  + (H20‘^)] 

(4.3e) 

+ 

II 

II 

NJ 

0 

+ 

II 

0 

(4.3f) 

" a-f)[(0H'^)  + (h'^)  + (H20'^)  ] 

(4.3g) 

where  the  parentheses  indicate  initial  concentrations,  and  brackets 
indicate  post-recombination  concentrations.  The  various  levels  of 
clustering  of  the  hydrated  proton  are  indicated  collectively  by 

By  post-recombination  is  meant  for  times  up  until  the  neutral 


reactions  get  underway. 
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When  we  think  of  Equations  (4.3a-g)  in  the  yield  rate  sense,  it  is 
necessary  to  include  the  yield  rate  for  H^O,  which  follow  from  mass 
conservation: 


(H^O)  = -(OH)  - (0)  - (OH^)  - (4.4a) 

= - [(H^)  + (H^O"^)]  -|[(H)  + (h"^)  + (OH)  + (OH"^)].  (4.4b) 

In  the  post-recombination  phase  additional  H2O  molecules  have  been 
broken  up  as  a result  of  electron  transfer  from  H2O  to  H'*’  and  OH"*" 
(reactions  b and  d)  to  yield  H2O'*’  which  then  loses  a proton  to  H2O 
(reaction  a)  to  yield  OH  and  H^O'*’.  In  all  the  other  reactions  H2O 
merely  serves  as  a catalyst  to  reactions  which  do  not  break  up  H2O 
molecules.  It  follows  that 

[HjO]  . (H^O)  - (Hh  - (OHh  (4.5) 


independent  of  f. 

We  may  check  our  predictions  for  the  post-recombination 
concentration  levels  against  the  results  of  Rio's  code,  both  for  the 
case  of  initial  pulse,  and  for  the  case  of  a steady  beam,  which 
integrates  out  from  the  initial  pulse  case — as  far  as  the  ion  reactions 
are  concerned.  We  find  that  our  predictions  are  borne  out  to  four 
significant  figures.  This  is  not  too  surprising,  since  all  we  were 
calculating  were  final  states,  or  asymptotic  limits. 
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From  the  above  it  is  clear  that  all  concentrations  of  interest 
except  for  the  individual  breakdown  of  the  collectively  defined  hydrated 
proton,  and  the  concentrations  of  the  unreacted  primary  ions,  which  are 
minuscule,  can  be  calculated  solely  on  the  basis  of  neutral  chemistry 
and  electron-proton  recombination. 

If  the  electron  yield  is  everywhere  greater  than  the  ion  yield, 
then  the  problem  reduces  to  one  of  only  neutral  chemistry.  If  this  is 
not  the  case,  then  the  mobility  of  the  electron  comes  into  play. 


Charge  Separation  and  Electron  Mobility 

The  charge  separation  which  is  present  already  in  the  initial 

deposition,  in  greater  or  lesser  degree,  depending  on  Incident  energy 

creates  an  electrostatic  force  field  which  Influences  the  motion  of  the 

electron  and,  to  a lesser  extent  because  of  their  far  greater  masses, 

the  motion  of  the  ions.  Specifically,  it  gives  rise  to  a drift 

velocity,  v,  given  by  v = pE  where  p is  defined  as  the  mobility 

coefficient  and  is  given  by  p = (e/m)v  where  v is  the  frequency  of 

c c 

collisions  (Allis,  1956;  Brown,  1958). 

A convective  term  is  thereby  Introduced  into  the  diffusion  kinetic 
equation  for  the  electron  (the  convective  term  for  the  positive  ions  can 
be  relatively  ignored); 


8c 
e 

8t 


pE  • Vc 

e 


2 

D V c 
e e 


kc  c 
P e 


(4.6) 
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where  Cp  is  the  concentration  of  positive  ions,  which  may  be  taken  to  be 
the  concentration  of  hydrated  protons  after  the  initial  yield  positive 
ion  species  have  reacted  with  H2O.  For  the  case  of  spherical  geometry 
Eq.  (4.6)  becomes 

+ JiM  ^ JL  z'  2 ^‘^e. 

9t  3r  9r  ^2  3r  9r  ^ ~ ^‘^p'^e  ('^•7) 

where  the  electrostatic  potential,  $(r),  is  given  by 

The  problem  of  ambipolar  diffusion  with  recombination  expressed  by  Eq, 
(4.6)  is  an  interesting  and  varied  one,  to  which  much  attention  has  been 
given  recently  (Christophorou,  1981).  A full  treatment  of  this  problem 
will  not  be  attempted  here.  In  its  place,  one  extreme  limit  of  the  two 
which  bracket  the  exact  solution  will  be  taken  as  a working 
approximation,  so  that  we  may  get  on  with  things.  However,  this  would 
be  an  interesting  future  direction  in  which  to  take  this  work. 

These  two  extreme  limits  consist  of  setting  p equal  to  0,  on  the 
one  hand,  and  on  the  other.  In  the  first  case,  the  electron  is  not 

influenced  by  the  field,  and  only  diffuses.  In  the  second  case,  the 
influence  of  the  field  is  overwhelming,  and  diffusion  is  bypassed  by 
instantaneous  global  recombination.  The  latter  course  is  chosen  here. 

We  therefore  set  f=l  in  Eqs.  (4.3a)  and  (4.3g),  and  invalidate  Eq. 

(4.3e). 

The  best  approximation  to  the  true  time  course  of  the  concentration 
of  the  thermal  electron  that  can  be  had  under  these  circumstances  would 


129 


probably  be  to  simply  renormalize  the  initial  yield  by  the  factor 
1/N(Eq),  and  allow  these  electrons  to  simply  diffuse. 

For  this  purpose  we  then  need  the  diffusion  coefficient  of  the 

thermalized  electron  through  H2O.  Rio  (Ph.  D.  thesis,  1983)  observed 

that  it  may  be  estimated  from  the  elastic  cross  section,  extrapolating 

low  energy  cross  section  data  of  Seng  (1975)  below  .35  eV. 

Specifically,  the  mean  free  path  between  collisions,  <£>,  is  given  by 

<0  = 1/a^^n.  The  diffusion  coefficient  is  given  exactly  by 
2 

= <v  /3v^>  = <v£/3>  (Brown,  1958)  where  is  the  collision 
frequency,  following  from  its  definition  as  the  proportionality  constant 
between  the  total  particle  flux  and  the  concentration  gradient.  The 
ensemble  average  can  be  approximated  by  neglecting  correlation  in  one  of 
two  ways  (Reif,  1965): 

« <v><£>/3  = (l/3)(8kgT/Tnn^)^/2^l/a^^n)  (4.9.1) 

(dimensional  estimate) 

« <v2><£>/3<v>  = (l/3)(3kgT/m^)(l/a^_^n)/(8kgT/7rm^)l/2 

(random  walk  estimate) 

These  two  estimates  differ  by  less  than  the  uncertainty  in  the 
extrapolated  estimate  of  A temperature  and  density  dependent 

value  for  D follows, 
e 

Knowing  D^,  we  obtain  p via  the  Einstein  relation, 

= (m/e)<v^/3>p 


(4.10) 
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(Allis,  1956;  Brown,  1958).  It  follows  then  that  ji/e  = D /kT.  The 

e 

general  problem  of  charge  separation  posed  by  Eq.  (4.7)  can  now  be 
addressed  quantitatively.  This  problem  is  left  for  future  work. 


Ion  Reaction  Time  Scales 

We  here  attempt  to  undertake  a quantitative  treatment  of  the  time 
scales  of  the  ion  reactions  and  their  dependence  on  medium  density  and 
injection  rate.  Alternatively,  we  could  simply  run  a kinetic  code,  such 
as  Rio's,  from  which  Fig.  4.1  was  obtained,  for  a few  selected  values  of 
medium  density  and  injection  rate,  and  find  out  this  way  without  any 
effort  for  either  the  initial  pulse  or  the  continuous  pulse  case.  We 
choose  to  go  the  first  route  out  of  interest,  since  the  latter  one  is 
always  available  to  us,  if  we  need  it. 

In  order  to  simplify  things,  we  may  first  of  all  ignore  the  primary 
species  OH'*’  and  which  are  produced  at  about  one-tenth  the  level  of 
H20''’.  We  may  therefore  Ignore  reactions  b-d.  Now,  only  considering 
reactions  a,  e,  f,  and  i,  g and  h,  we  may  further  approximate  by 
combining  the  various  hydrated  proton  species  into  one  average 
species.  We  are  left  then  with  reaction  a with  rate  constant  unchanged, 
and  an  average  of  reactions  i,  g and  h to  which  we  assign  an  average 
rate  constant.  At  lower  medium  densities,  this  rate  should  be  shifted 
down  toward  the  the  singly  hydrated  case  (the  hydroniura  ion).  At  higher 
medium  densities  this  rate  should  be  shifted  up  toward  the  triply 
hydrated  proton. 

The  proton  transfer  reaction  (a)  can  be  treated  as  a single  body 
decay  because  of  the  constant  background  of  H2O.  The  time 


constant  is 
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found  via  Eq.  3.6  in  the  case  of  initial  pulse  or  Eq.  3.3  in  the  case  of 

continuous  pulse.  For  the  system  of  Fig.  3.1  (initial  pulse)  this  turns 

*“8 

out  to  be  about  10  sec,  which  agrees  with  the  figure.  (Although  not 
shown,  = 10  mole/cc,  as  can  be  seen  by  noting  the  time  course  of 

H^O  and  the  rate  constant  for  reaction  a.) 

The  reaction  time  of  the  electron— proton  recombination  is  found 
either  from  Eq.  (3.9)  in  the  initial  pulse  case  or  from  Eq.  (3.22)  in 
the  continuous  pulse  case.  Equation  (3.9)  (or  its  asymptotic  case,  Eq. 
(3.6))  does  not  spell  out  reaction  time.  We  may  do  that  here,  using  the 
e-folding  interval  as  a definition.  One  obtains 


T 

C 


= - <=20> 

■< 

^10  ~ ^*^20 

^^20  ^ ‘^lO  ~ ^^^20 

(4.11) 


For  the  system  of  Fig.  4.1  the  average  rate  implies  an  e-folding  time  of 
“8 

6 X 10  sec.  This  is  somewhat  low.  Using  the  rate  for  instead, 

considering  the  low  density,  one  obtains  1.3  x 10~^  instead,  which  is 
in  better  agreement. 

The  above  analysis  applies  to  the  case  in  which  the  reaction  time 
of  recombination  is  later  than  the  reaction  time  of  proton  transfer.  In 
this  case  the  reaction  time  of  recombination,  the  last  step,  defines  a 
close-out  time.  In  the  opposite  case  the  proton  transfer  reaction  is 
the  rate  limiting  step',  and  defines  both  onset  and  close  out. 

It  is  convenient  to  leave  this  question  of  the  time  scales  of  the 
ion  reactions  momentarily  in  order  to  go  on  to  a discussion  of  the 
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neutral  reactions  in  the  next  chapter.  However,  we  wili  return  to  it  at 
the  end  of  that  chapter  and  complete  our  discussion  of  it  there. 


CHAPTER  V 

KINETICS  OF  FULL  NEUTRAL  REACTION  SET 
Introductlog 

The  radiological  context  dictates  four  degrees  of  freedom:  medium 

density,  injection  rate,  incident  energy  and  temperature.  Room 
temperature  is  the  most  natural  choice,  especially  for  the  liquid  state. 
In  this  connection  it  is  important  that  we  make  it  clear  that 
calculations  for  liquid  density  based  on  spatial  yields  computed  from 
gaseous  phase  cross  sections  and  on  a gas  phase  reaction  set  are  only  to 
be  considered  as  extrapolations. 

Taking  as  our  point  of  comparison  the  spatial  yields  at  r=0  as  we 
had  done  in  Chapter  III,  two  of  these  degrees  of  freedom,  namely 
incident  energy  and  injection  rate,  as  a result  of  the  near  constancy  of 
g-values  for  E^  > 100  eV,  combine  effectively  into  one  degree  of 
freedom  for  this  corresponding  kinetic  system— the  "strength",  or 
scaling,  of  the  source  terras  of  the  assorted  species. 

If  it  is  a true  result  that  the  range  of  the  high  spectral  plume 
relaxes  as  E^  increases,  then  this  fact,  along  with  the  clear  decrease 
in  the  high  spectral  plume  g-value  with  increasing  E^,  can  be  taken  to 
imply  essentially  single  range,  or  single  exponential,  spatial  profiles 
of  the  yields  at  higher  E^,  say,  E^  > 1 keV.  The  general  picture  of 

spatial  profiles  of  the  yields  which  emerges  is  then  one  based  only  on 
the  total  g-value. 
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For  the  post-recombination  yields,  we  obtain,  from  Fig.  2.20,  the 
following  relative  values  of  the  total  g-values: 

% ” ®0H’  ®0  ” 

The  temperature  of  400  K is  somewhat  arbitrarily  chosen  to  work  at. 
For  the  liquid  density  case  both  400  K and  300  K are  calculated  at. 
Because  of  the  energetics  of  the  reactions  in  this  temperature  region, 
this  does  not  lead  to  a significant  difference  in  the  character  of 
results.  That  this  should  be  so  can  be  already  appreciated  from 
considering  Table  4.2  where  reaction  rates  at  300  K as  well  as  400  K are 
presented  (the  general  temperature  dependences  are  given  in  Westbrook, 
1981).  This  can  be  seen  for  the  third  body  reactions  as  well  as  for  the 
two  body  reactions. 

This  reaction  set.  Intended  to  apply  to  combustion,  is  more 
extensive  than  we  really  need  for  the  radiological  problem  in  the 
vicinity  of  room  temperature.  Most  of  the  reactions  can  be  taken  as  one 
way  in  this  temperature  domain.  It  is  easy  to  see  that  this  is  so  for 
all  the  third  body  reactions.  That  it  is  so  for  most  of  the  two  body 
reactions  can  be  seen  from  Table  5.1  where  the  ratio  of  forward  rate 
coefficient  to  backward  rate  coefficient  or  the  reverse  is  given — (the 
ratio  of  the  larger  to  the  smaller  is  given).  Reaction  2 has  the  most 
equal  energetic  balance.  The  reactions  for  which  the  order  of  magnitude 
of  the  ratio  is  below  10  may  be  driven  in  the  reverse  direction  if  the 
concentrations  are  just  right.  This  however  is  the  exception.  Those 
for  which  the  order  of  magnitude  is  higher  yet  can  be  taken  as 
effectively  one  way  reactions.  A minimal  reaction  set  which  is  a subset 


135 


TABLE  5.1 

Ratios  Between  Forward  and  Backward  Rate  Coefficients 
for  the  Two-Way  Reactions 

(Ratios  are  of  the  larger  rate  coefficient  to  the  smaller.) 


Reaction 

Ratio  at  300  K 

at  400  K 

1 

4.3(10) 

5.3(7) 

2 

12.0 

5.3 

3 

3.5(11) 

2.5(8) 

4 

2.5(10) 

4.4(7) 

5 

1.3(22) 

3.0(16) 

8 

2.5(40) 

1.9(30) 

9 

1.4(29) 

1.6(22) 

10 

1.8(41) 

7.3(30) 

11 

9.6(51) 

5.2(38) 

12 

5.5(29) 

1.4(22) 

14 

1.8(11) 

3.5(8) 
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of  this  one  can  be  identified.  Rather  than  bother  with  this  at  this 
point,  we  will  let  the  kinetic  calculations  tell  us  which  ones  are 
contributing,  since  little  extra  expense  is  involved  in  the  numerical 

calculation  from  more  reactions,  this  coming  primarily  from  number  of 
species. 

It  should  be  noted  that  the  favored  one-way  reactions  are  all 
either  two  body  or  effectively  two  body  (third  body  is  primarily 
background  H2O).  The  two  effectively  one  body  decay  reactions  (3f  and 
4f)  are  energetically  unfavored. 


Principal  Pathways  and  Reaction  Times 

Setting  all  yields  (i.  e.  source  terms )~f or  H,  OH,  and  0— to  be 

of  the  form,  y^  = = (3^^d^IZ^/87t)e  ^ » where  I is  the 

Injection  rate,  since  we  want  to  calculate  pure  kinetics  for  the  yields 

obtaining  at  r=0,  we  turn  our  attention  on  rather  than  Z^,  the 

spatially  integrated  yield.  We  therefore  define  a quantity,  S = Id^ 

which  we  may  call  the  source  strength. 

We  are  interested  in  the  variation  in  the  kinetic  time  courses  with 

medium  density,  d,  and  source  strength,  S.  We  define  dimensionless 

density,  injection  rate  and  source  strength  in  terms  of  convenient 

reference  values:  h = d/d  , i = I/l  and  g = s/S  where  S = I d ^ 

^ o 000 

where  = 10^2  ///sec  and  d^  = 10~^  gm/cc  = 6.02xlol5  #/cc. 

We  now  investigate  the  g-h  grid.  The  grid  points  (g,h)  calculated 
for  are  labelled  in  Fig.  5.1.  The  kinetic  time  courses  for  those 
labelled  by  an  alphabetic  character  other  than  an  o are  presented  in 
Figs.  5.2a-i,  Fig.  5.3j-k,  Fig.  5.4£  and  Fig.  5.5m.  These  are  for  post- 
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Figure  5.1.  Medium  density-source  strength  grid.  Overlaid  with 
isocontours  of  (1)  injection  rate  for  equivalent 
yield  at  r = 0,  (2)  kinetic  onset  time,  (3)  borderline 
of  fast/slow  chemistry,  and  (4)  borderline  of  depletion 
regime  for  diffusion  only.  Points  at  which  kinetic  runs 
were  made  are  labelled. 


Figure  5.2.  Kinetic  time  course  of  post-recombination  neutral 
system  for  T = 400  K,  = 1 keV,  yield  rates 
obtaining  at  r = 0,  for 

(a)  g = lO,,  h = 1 (t  = .015  sec), 

(b)  g = 10^,  h = 10  ?T  = .0015  sec), 

(c)  g = 10^,  h = 100  (t  = .0015  sec), 

(d)  g = 1.4  X 10^.  h = 10^  ^ (t  = .0015  sec), 

(e)  g = 1.0  X 10j\  h = 3067  (t^  = 1.3  x 10  sec) 

(f)  g = 1.4  X 10^,  h = 10^  (t  =*^.0015  sec), 

(g)  g = 1.4  X 10^,  h = 10^  = .0015  sec), 

(h)  g = 1.4  X 10^,  h = 10  (t^  = .00015  sec), 

(i)  g = 2.5  X lO-*,  h = 5.57  x^lO  (x  = .00015  sec) 
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Figure  5.3.  Kinetic  time  courses  below  the  Figure  5.3.  Kinetic  time  courses  below  the 

minimum  g for  convergence  of  pattern:  minimum  g for  convergence  of  pattern 

(j)  g = 140,  h = 10  (k)  g = 1.4,  h = 10^ 

(t,-  = .0015  sec).  (x  = .015  sec). 
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Figure  5.4.  Kinetic  time  course  below  the 
minimum  g for  convergence  of 
pattern: 

U)  g = 250,  h = 5.57  X 10® 

(t  = .0015  sec). 
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Figure  5 


.5.  ^netic  time  course  illustrating  depletion: 
(m)  g = 10°.  h = 1 (x^  = .0015  sec). 
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recombination  spatial  yields  at  T = 400  K and  = 1 keV.  Fig.  5.6 
presents  the  same  system  as  in  Fig.  5.21  (liquid  density)  except  with 
the  difference  that  it  is  at  300  K.  Systems  (a)-(i)  with  the  exception 
of  (e)  all  have  = .015  sec,  .0015  sec  or  .00015  sec.  The  choice  of  g 
in  system  (e)  is  made  to  agree  with  the  diffusion  kinetic  system  of  H + 
OH  in  Chapter  III  and  with  the  second  injection  rate  for  the  second 
density  considered  in  Chapter  VI.  System  m Illustrates  depletion. 

Choices  (j)-(k)  and  (1)  Illustrate  a significant  effect  which  we  did  not 
Investigate  further.  This  is  the  following.  In  all  other  cases  we 
observed  no  variation  of  pattern  of  the  kinetic  time  courses  with  g— 
only  a rescaling  of  concentration  and  time  scale.  However  at  these  low 
source  strengths,  we  observe  the  pattern  begin  to  change.  This  would  be 
an  Interesting  avenue  of  further  study. 

In  the  g-h  grid  in  Fig.  5.1  isocontours  of  i are  plotted  from  the 
equation,  g = hi  . Thus,  noting  maximum  and  minimum  physically 
meaningful  values  of  i,  we  limit  the  region  of  the  g-h  grid  which  we  are 
interested  in.  The  mathematical  properties  of  grid  points  with  respect 
to  the  associated  kinetic  systems  are,  however,  independent  of  these 
physical  considerations.  The  meaning  of  unchanging  pattern  with  respect 
to  g above  a minimum  value  can  be  clarified.  This  means  that  only  the 
vertical  and  time  axes  need  by  rescaled* 

The  evolution  of  pattern  seen  in  Figs.  5.2a-l  is  now  discussed.  We 
first  look  at  the  early  time  regime.  Here  chemistry  is  still  only 
feedforward  to  the  product  species,  and  the  reactions  which  produce  them 
are  noted:  lb  for  O2,  13b  for  H2O2,  and  primarily  7f  for  HO2.  We  are 
more  interested  in  the  first  two.  Since  reaction  13b  is  third 


body,  the 
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Figure  5.6. 


Kinetic  time  course  for  same  yield  strength  and 
density  (liquid)  as  in  system  (i)  but  at  T = 300  K 
(system  (n)).  Kinetic  time  course  for  system  (i)  is 
overlaid  (dashed  curve). 
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production  rate  of  H2O2  increases  with  increasing  h.  H2O2  production, 
below  that  of  O2  at  low  densities,  is  considerably  above  it  at  the 
higher  densities. 

It  is  noteworthy  that  H is  fairly  unreactive  at  low  densities.  The 
two  body  reactions  it  is  involved  in  are  not  very  fast.  At  higher 
density,  three  body  (radical  recombination)  reactions  are  favored,  and  H 
is  more  reactive  in  this  density  regime. 

We  now  home  in  on  three  representative  cases:  system  (a),  the  low 
density  gas;  system  (e),  the  1 atm  gas;  and  system  (i),  the  liquid 
density  case,  essentially  the  extremes  and  about  the  midpoint.  At  each 
of  these  densities,  we  plot  the  net  production  rate  time  courses 
(including  source  terra).  According  to  the  kinetic  equations,  the  net 
production  rate  should  equal  the  time  derivative  of  the  species 
concentration.  This  will  be  checked.  These  plots  are  given  in  Fig.  5.7 
for  system  (a).  Fig.  5.9  for  system  (e),  and  Fig.  5.10  for  system  (1). 

In  order  to  identify  the  dominant  reactions,  also  known  as  the 
principal  pathways,  we  have  also  plotted  time  courses  of  Individual 
reaction  rates  (i.  e.  the  net  two  body  or  the  associated  net  effective 
two  body  reaction  rate).  These  are  given  in  Figs  5.8,  5.11  and  5.12  for 
systems  (a),  (e),  and  (i),  respectively. 

Finally,  these  reaction  rates  are  multiplied  by  the  appropriate 
signed  stoichiometric  factors  to  give  production  rates  of  species  per 
reaction.  Time  courses  of  these  quantities  are  plotted  in  Figs.  5.13, 
5.14  and  5.15  respectively  for  systems  (a),  (e),  and  (1).  The  assignment 
of  subfigures  is  the  same  in  all  of  these,  namely:  (a)  H,  (b)  OH,  (c) 

^2*  (d)  0,  (e)  HO2,  (f)  O2,  (g)  H2O2  and  (h)  H2O. 
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Figure  5.11.  Time  course  of  reaction  rates  for  kinetic  system  (e) . 
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Figure  5.12.  Time  course  of  reaction  rates  for  kinetic  system  (1) . 


Figure  5.13.  Time  courses  of  production  rates  of  (a)  H,  (b)  OH, 
(c)  H2,  (d)  0,  (e)  HO2,  (f)  O,,  (g)  H2O2,  and 
(h)  H2O  for  kinetic  system  (a)  from  each  reaction. 
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Figure  5.13  (h) 


Figure  5.14.  Time  courses  of  production  rates  of  (a)  H,  (b)  OH, 
(c)  H2,  (d)  0,  (e)  HO2,  (f)  O2,  (g)  H2O2,  and 
(h)  for  kinetic  system  (e)  from  each  reaction. 
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Figure  5.15.  Time  courses  of  production  rates  of  (a)  H,  (b)  OH 
(c)  H2,  (d)  0,  (e)  HO2,  (f)  O2,  (g)  H2O2,  and 
(h)  H2O  for  kinetic  system  (i)  from  each  reaction 
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From  this  series  of  time  courses  we  can  precisely  identify  the 
principal  pathways  overall  and  in  the  production-destruction  of  each 
species.  Let  us  now  describe  the  more  significant  features.  We  may 
start  with  the  evolution  of  reaction  6,  H + OH  recombination,  and  its 
weaker  companion,  reaction  13,  OH  + OH  recombination.  In  system  (a), 
reaction  13  is  barely  noticeable,  while  6 is  of  average  importance  among 
the  principal  pathways  at  the  later  stage  (following  the  second  reaction 
time  in  this  system).  In  systems  (e)  and  (i)  they  are  both  the  strongest 
two  reactions  in  system  (1)  for  all  time  seen,  in  system  (e)  until 
later  time  when  a second  group  of  reactions  start  up  (similar  to  the 
situation  in  system  (a)).  In  system  (a),  two  reactions  which  are 
relatively  weak  in  (e)  and  (i)  pre-stage  the  show  by  two  orders  of 
magnitude  in  time.  These  are,  namely,  reactions  1 and  3. 

We  may  now  calculate  reaction  time.  First  we  look  at  low  density. 
The  reactions  1 and  3 form  a feedback  loop  on  OH.  The  rate  coefficient 
of  reaction  3 is  lower  by  a factor  of  10  than  that  of  reaction  1. 

However  the  ratio  of  source  terms  y^^^/y^  = 8 mostly  makes  up  for  it.  As 
a result  OH  scavenges  itself  in  the  net  balance  about  as  fast  as  it 
scavenges  the  less  plentiful  0,  and  OH  is  seen  to  fall  off  in  sync  with 
0,  and  then  continues  after  the  source  term  supply  of  0 has  already  been 
neutralized.  The  reaction  time  is  given  in  Eq.  (3.22)  for  this 
system.  We  then  have 


X 

c 


= '/^"ib^^H  - ^0>  • 


(5.1) 


At  higher  density,  the  reaction  time  is  that  of  reaction  6,  and,  again 
following  Eq.  (3.22),  we  obtain 
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T 

C 


(5.2) 


Substituting  in  numbers,  this  comes  out  to  t = l.5/g^^^  and 

1 /2 

Tc  = 24. /(gh)  . These  reaction  times  intersect  at  h = 256.  The  h 

value  at  which  early  time  H2O2  and  O2  production  are  equal  is  140. 
Contours  of  constant  are  plotted  in  Fig.  5.1. 

In  the  diffusion  kinetic  case,  it  is  necessary  to  know  also  the 
diffusion  time,  r^.  Recalling  from  Chapter  III  that 

Td  = = l/Wd^  h (D  = D^(d^/d),  we  note  that  for 

H,  Dq  = 6400  cm^/sec,  for  OH,  = 1830  cm^/sec  (Table  6.1).  We  may 
take  an  average  of  4000  cm^/sec  as  a fuzzy  value  in  the  H-OH 
recombination  region.  In  the  low  density  region,  noting  that 

°0  ” °0H*  ” °0-  ®^o  = 1 6d^  = .25  cm~^ 

therefore  = .004/h  sec  in  the  first  case  and  .008/h  sec  in  the 
second.  The  dependence  of  the  ratio  x^/x^  on  g and  h follows  at 
once.  The  isocontour  = 1 is  plotted  in  Fig.  5.1. 

We  may  also  consider  depletion.  The  time  at  which 

~ ~ ^dpl  found  to  be  equal  to  10^  h/g  (substituting  in 

numbers  from  the  G values  given  in  Chapter  II).  Comparing  against 
diffusion  convergence  time  in  the  low  density  case  where  diffusion 
depletion  tends  to  outpace  chemistry,  this  time  is  significant. 

However,  in  the  high  density  region,  this  time  is  no  longer  Important 
because  reaction  6 quickly  returns  much  or  most  of  radiolyzed  H2O  to  its 
original  state.  The  isocontour  = 1 Is  also  plotted  in  Fig. 

5.1. 

The  isocontours  plotted  in  Fig.  5.1,  as  has  been  noted,  are  for 
Eq  = 1 keV  and  T = 400  K.  The  complete  picture  of  and  T dependence 
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can  be  seen  from  a representative  set  of  these  g-h  grids,  overlaid  with 
isocontours,  over  the  - T grid. 


Particle  Conservation  and  the  Stiffness  Problem 

Kinetic  systems  of  equations  are  famous  for  the  stiffness  problem 
which  results  if  disparate  time  scales  are  involved  (Dahlquist,  Bjorck, 
1974;  Burden,  Faires  and  Reynolds,  1978).  That  can  easily  be  seen  to  be 
the  case  here  in  systems  (a)  and  (e)  a second  reaction  time  manifests 
itself,  while  in  systems  (e)  and  (i)  reaction  rates  are  seen  to  divide 
into  two  groups — a very  dominant  few  and  a group  of  several  small 
ones.  In  system  (e),  this  group  increases  in  importance  at  later  times. 

In  order  to  check  the  accuracy  of  our  calculations  we  have  done  two 
things.  Certainly,  this  is  only  a beginning.  One  was  to  plot  the  net 
production  rates,  so  that  we  can  compare  them  with  the  slopes  of  the 
concentration  time  courses.  The  other  is  that  we  calculated  the 
quantity  (H  atom)/2(0  atom)  - 1,  which  measures  mass  balance  (but  not 
necessarily  mass  conservation).  This  quantity  is  graphed  in  Fig.  5.16 
for  kinetic  systems  (a),  (b),  (c),  (f),  (h)  and  (i).  Contribution  from 
H2O  species  was  not  taken  into  account.  In  the  high  density  cases, 
since  so  much  of  the  initial  yield  of  H and  OH  is  recombined,  this 
exaggerated  the  reading. 

In  all  cases,  either  because  of  overproduction  of  H atoms  or 
underproduction  of  0 atoms  (more  likely  the  case  because  of  the 
levelling  off  of  the  reactive  radicals  OH  and  0) , this  ratio  increases 
Bft&v  reaction  time.  It  peaks  at  12—15%  for  the  higher  densities,  and 


Figure  5.16.  Time  courses  of  mass  conservation  in  numerical 

calculation.  (H  atom)/(0  atom)  - 2 is  the  quantity 
plotted,  where  H2O  species  is  not  counted.  For 
kinetic  systems  ta)»  (b),  (c),  (f),  (h)  and 
(i)  in  that  order. 
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Figure  5.16  (a) 


Figure  5.16  (b) 
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Figure  5.16  (d) 
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Figure  5.16  (h) 


Figure  5.16  (i) 
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rises  like  a singularity  at  the  low  densities,  but  at  the  very  later 

time  after  the  (first)  reaction  time  it  starts  its  rise  slowly. 

Comparing  the  plots  of  net  production  rates  with  the  concentration 

time  courses,  we  observe  discrepancies,  mostly  at  later  times.  These 

later  times  are  it  should  be  noted rather  ambitious.  Rio  (Ph.  D. 

thesis,  1983;  see  Fig  4.1,  this  thesis)  ran  his  kinetic  code  for  the 

full  ion-reaction  system  for  initial  pulse  source  term  out  to  about 

~ 200-1000.  We  are  running  here  for  steady  source  term  out  to 

t/T  = 40,000. 
c 

It  is  of  interest  to  note  here  (Rio,  private  communication,  1983; 

follow  up  by  author)  that  the  kinetic  code  he  used,  which  he  had  adapted 

from  a previous  code  developed  by  Walker,  Pamidimukkala,  Meyreddy  and 

Green,  which  is  based  on  the  Gear  method  (Hindmarsh,  1974)  begins  to  act 

up,  when  calculating  for  steady  source  term,  around  t ==■  t . The 

c 

results  it  give  are  strange:  the  concentrations  of  the  product  species 

start  to  become  jagged  in  their  time  courses  and  even  assume  negative 
velues  intermittently  which  in  their  absolute  value  are  nearly 
continuous  with  their  general  trends. 

Returning  to  the  plots  of  net  production  rates,  we  consider  two 
examples,  the  HO2-O2  subsystem  in  system  (i)  (it  is  incipient  in  system 
(e))  and  the  time  course  of  HO2  in  the  low  density  case.  Looking  to  Fig. 
5.10  in  the  first  example,  we  see  the  net  production  rate  of  O2  go 
negative  from  initial  value  of  zero  (while  that  for  HO2  does  the 
opposite).  This  conflicts  with  the  positive  definiteness  of  its 
concentration  time  course.  From  Fig.  5.15e  and  f it  is  apparent  that 
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reaction  7 is  the  source  of  the  problem.  We  also  notice  the  unlikely 
way  in  which  reaction  rate  7 cuts  across  the  other  reaction  rates  in 
Fig.  5.12. 

Looking  at  Figs.  5.7  and  5.13e  in  the  second  example,  we  see  the 
net  production  rate  of  HO2  fall  from  an  initial  value  of  zero, 
conflicting  with  the  positive  definiteness  of  the  concentration  time 
course.  We  see  the  problem  in  Fig.  5.13e — the  net  production  rate  is  the 
small  difference  of  large  quantities.  Both  of  the  above  examples  are 
examples  of  stiffness. 

The  consolations  are  (1)  that  the  discrepancies  occur  for  the 
species  in  lower  concentration,  and  (2)  that  the  overall  trends  or 
dominant  trends  are  calculated  to  an  adequate  level  of  accuracy,  and  (3) 
the  calculation  is  numerically  stable  so  at  least  some  answer  is 
given.  That  this  is  so  is  a result  of  the  implicit,  linearized 
calculation  of  the  nonlinear  kinetic  production-destruction  terms  (see 
Appendix  A),  which  is  known  to  be  numerically  stable  (Dash  and 
Pergament,  1978;  Dahlquist  and  Bjorck,  1974). 

There  is  perhaps  a fourth  consolation.  Our  reason  for  interest  in 
the  kinetics  is  to  base  a comparison  with  diffusion  kinetics.  The 
important  question  then  becomes  — how  does  diffusion  affect  the 
stiffness  problem?  Since  the  effect  of  diffusion  is  to  slow  down  the 
chemistry  by  placing  an  upper  bound  on  the  growth  of  the  concentration 
levels  of  the  reactants,  it  would  seem  reasonable  to  conjecture  that  the 
effect  of  diffusion  on  the  stiffness  of  the  kinetic  system  is  one  of 
roftig^tion.  It  would  be  interesting  to  attack  this  question  rigorously. 
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Temperature  Sensitivity 

In  order  to  exhibit  the  effect  of  dependence  on  temperature,  the 
kinetic  time  course  of  the  same  kinetic  system  as  system  (i)  (liquid 
density  at  400  K)  only  with  temperature  changed  to  300  K (this  system 
will  be  denoted  as  system  (n))  has  been  plotted  in  Fig.  5.6.  The 
kinetic  time  course  of  system  (i)  is  overlaid  for  the  sake  of  easy 
comparison.  The  sensible  differences  between  the  two  systems  are 
visible  from  this  figure.  Overall,  taking  into  account  the  distortion 
of  log  vertical  axis,  the  general  trends  are  unaffected. 

Specifically,  we  note:  (1)  initial  H2O2  production  increased  by  a 
factor  of  about  eight  (reflecting  the  increase  in  rate  constant — see 
Table  4.2),  (2)  somewhat  earlier  reaction  time,  (3)  lower  final  levels 
of  H,  0,  and  OH,  decreased  by  factors  of  2,  0.5,  and  2.5,  respectively, 
^^itial  ©2  and  HO2  production  rates  are  unaffected.  Later  time  ©2 
production  rate  is  also  unaffected.  However,  later  time  HO2  ambient 
level  is  down  by  factor  of  2. 

Net  production  rates  are  graphed  in  Fig.  5.17,  reaction  rates  in 
Fig.  5.18,  and  reaction-specific  production  rates  in  Fig.  5.19,  in  like 
manner  as  the  three  systems  previously  detailed.  Two  distinct  energetic 
effects  are  noted:  (1)  most  third  body  radical  recombination  reactions 
are  speeded  up  (two  are  unchanged),  and  (2)  all  two  body  reactions  are 
slowed  down  in  both  directions,  however  with  the  preferential  sense 
retained,  and  the  general  order  of  magnitude  neighborhood  of  the 
preferential  sense  also  retained  (see  Table  5.1). 

The  die  off  of  two  body  reactions  is  most  strikingly  illustrated  by 
the  cases  of  reactions  3 and  4,  which  have  completely  disappeared.  The 
increases  in  third  body  reaction  rates  are  far  from  proportional.  This 
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Figure  5.17.  Time  coure  of  net  production 
rates  for  kinetic  system  (n) . 
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Figure  5.18.  Time  course  of  reaction  rates  for  kinetic  system  (n) . 


Figure  5.19.  Time  course  of  reaction-specific  production  rates 
of  (a)  H,  (b)  OH,  (c)  H2,  (d)  0,  (e)  HO2,  (f)  O2 
(g)  H2O2,  and  (h)  H2O  for  kinetic  system  (n). 
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is  most  notable  in  the  case  of  H + OH  and  OH  + OH  (reactions  6 and 
13).  The  ratio  goes  from  0.6  at  400  K to  3.0  at  300  K.  The  effect  is 
significant:  reaction  time  now  is  keyed  to  the  latter  reaction  rather 

than  the  former,  and  the  transient  spurt  of  H2O2  production  is  doubled 
in  peak  height,  as  a result  of  which  initial  OH  destruction  comes 
primarily  from  H2O2  production,  whereas  at  400  K,  it  comes  equally  from 
H2O2  and  H2O  production.  This  can  be  seen  by  comparing  Figs.  5.15  (b) 
and  5.19  (b). 

Separation  of  Ion  and  Neutral  Time  Scales 
Now  let  us  pick  up  where  we  left  off  with  the  ion  reaction  time 
scales.  We  will  first  compare  time  scales  for  incident  energy  of  1 keV 
and  T = 400  K.  We  will  then  also  consider  T = 300  K.  Quantitative 
results  will  be  expressed  in  terms  of  the  convenient  units  of  source 
strength,  medium  density  and  injection  rate,  g,  h and  i,  defined  two 
sections  back.  By  use  of  the  quantity,  "source  strength,"  we  imply 
conditions  at  r=0. 

We  first  note  the  reaction  time  for  reaction  a (proton  transfer: 
H20^  + H2O),  given  by 

Tpt  = 8.8  X 10  ^/h.  (5.3) 


Next  we  note  the  reaction  times  for  reactions  i,  g and  h (hydrated 
proton-electron  recombination),  given  by 
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^ .0052/g^/^ 

(H30^) 

"e-p  = " 

.0036/g^/^ 

(h'^(H20)2) 

(5.4) 

L .0028/g^^^ 

(H^(H20)3) 

We  may  ask,  given  a particular  density,  for  what  physically 
possible  injection  rate  (i.  e.,  one  which  does  not  cause  local  depletion 
the  medium,  much  less  give  rise  to  local  heating  or  arcing)  does  the 
recombination  reaction  outpace  proton  transfer?  (This  question  was 
asked  before  in  Chapter  IV.  We  are  now  going  to  answer  it.)  To  answer 
this  question,  we  set  implementing  an  average  value  of 

T^_p,  which  implies 

9 2 

g=1.5xl0h,  equivalently. 


Because  of  the  geometric  factor,  this  injection  rate  decreases  with 
increasing  density.  For  the  case  of  liquid  density,  Eq.  (5.5)  implies 
an  injection  rate  of  2.7x10^'^  #/s  or  a current  of  17  ya.  The  reaction 
time  is  then  1.6x10  sec  for  both  reactions. 

reaction  time.  The  reaction  times  of  the  neutrals  at  400  K 
in  terms  of  g and  h are  given  in  section  1.  We  may  repeat  them  here: 


" 1.5/g^/^,  h < 256 
^ 24/(gh)^/^,  h > 256,  f = 1 . 


(5.6) 
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It  should  be  noted  that  the  reaction  time  for  h < 256  is  independent  of 
f,  the  degree  of  charge  recombination,  since  H radical  is  not 
involved.  The  reaction  time  for  h > 256  does  depend  on  f somewhat  at 
400  K,  but  not  at  300  K where  OH  + OH  recombination  outpaces  H + OH 
recombination.  At  400  K,  the  difference  between  the  H + OH  reaction  time 
for  f=l  and  for  f=0  can  be  detailed.  Referring  to  the  graphs  of  initial 
and  post-recombination  g-values  at  the  end  of  Chapter  II,  and  to  Eq. 
(4.3a),  we  find  that  the  yield  rate  for  H is  increased  by  a factor  of  3 
in  changing  f from  0 to  1.  As  a result,  we  find  that  the  reaction  time 
of  H + OH  is  decreased  by  a factor  of  l/Zl  (see  Eq.  (3.22)).  Since  OH 
would  move  in  to  take  the  place  of  the  missing  H in  the  case  f=0,  it 
follows,  noting  the  different  reaction  rate,  that  the  reaction  time  is 
delayed  only  by  20%,  not  40%.  At  300  K,  reaction  time  for  h > 256  is 
given  by 

- U/(gh)*^2_ 

For  liquid  density,  this  becomes 

- .0047/g‘/2  . 

Thus,  at  liquid  density  the  separation  of  ion  and  neutral  time  scales 
has  disappeared.  In  fact  the  close-out  time  of  charge  recombination 
exactly  dovetails  with  the  onset  time  of  radical  recombination,  that  is, 
for  i < 270,  equivalently,  I < 17  pa,  or  equivalently,  g < 4.6x10^^. 
Referring  to  Fig  5.1,  we  note  that  this  is  already  extremely  into  the 
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early  chemistry  region  of  the  g-h  grid:  the  early/late  chemistry  border 

stands  around  g = 10^^,  or  equivalently,  i = 6x10“^,  for  the  liquid 
density  (to  within  20%  for  T = 300-400  K) . 

From  the  above  it  is  clear  that  at  all  injection  rates  in  the 
liquid  case  adequate  time  separation  of  charge  recombination  and  the 
neutral  reactions  cannot  be  assumed.  At  much  higher  injection  rates, 
adequate  separation  in  time  of  proton  transfer  and  charge  recombination 
cannot  be  assumed,  either.  This  latter  point  is  not  so  interesting, 
since  such  injection  rates  are  not  typical.  The  previous  point, 
however,  provides  still  more  motivation  to  upgrade  the  code,  DIFKIN,  to 
account  for  local  charge  recombination~a  point  that  was  mentioned  in 
Chapter  IV. 

In  the  gaseous  case,  up  to  high  pressures  near  liquid  density,  the 
separation  in  time  scales  between  ion  and  neutral  reactions  is  there. 

In  this  case,  the  motivation  to  include  charge  recombination  arises  only 
from  spatial  charge  separation. 


CHAPTER  VI 

DIFFUSION  KINETICS  OF  FULL  NEUTRAL  REACTION  SET 
Representative  Cases 

Diffusion  kinetics  are  calculated  for  the  three  representative 

densities  chosen  in  Chapter  V,  at  one  of  which  diffusion  kinetics  of  the 

single  reaction  H + OH  + M — > H2O  + M was  considered  in  Chapter  III. 

As  in  effect  in  Chapter  III,  we  consider  three  injection  rates  at  each 

density:  for  which  we  have  fast,  borderline  and  slow  chemistry, 

specifically  = x^/10,  = x^,  and  x^  = lOx^.  This  is  with  the 

exception  of  the  low  density  at  which  only  two  cases  are  considered: 

slow:  X = lOOOx  and  somewhat  fast:  x = x /2.5. 

c d c d 

The  diffusion  coefficients  have  been  calculated  according  to 

Chapman-Enskog  theory  (Hirshf elder,  Curtiss  and  Bird,  1953;  Monchick  and 

Mason,  1961)  utilizing  the  molecular  parameters  compiled  by  Warnantz, 

Miller  and  Kee  (1983).  They  are  given  in  Table  6.1. 

A schematic  of  cases  (i,h,T),  equivalently,  (g,h,T)  and  a guide  to 

figures  is  given  in  Table  6.2.  These  cases  are  also  shown  as  points  on 

^ g~h  grid  in  Fig.  6.2,  fitted  out  with  isocontours  of  x , x /x  and 

c c d 

^dpl 

Before  comparing  these  systems,  it  is  necessary  that  the  first 

system  be  mapped  into  the  corresponding  system  for  which  x = lOx 

c d' 

Since  we  are  interested  only  in  the  relative  values  of  the 
concentrations  we  simply  multiply  the  product  species  concentrations  by 
the  appropriate  factors  that  bring  them  into  the  correct  relationship 
with  the  initial  yield  species.  Those  factors  are  1400  for  O2  and  H2O2 
and  2x106  foj-  HO2. 


192 


TABLE  6.1 

Binary  Diffusion  Coefficients  in  H2O 
(in  units  of  cm^/sec) 


species 

Do 

H 

1.292 

6395 

«2 

.8644 

4279 

0 

.3754 

1858 

OH 

.3695 

1829 

H2O 

.2502 

1238 

HO2 

.2339 

1158 

H,0, 

.2327 

1152, 

O2 

.2305 

1141, 

is  diff.  coef.  at  p = = 1 atm  and  T = T,  = 300  K 

(n  = n^  = 2.45x10^^  #/cc) 

Dq  is  diff.  coef.  at  n = 6x10^^  #/cc  and  T = 400  K,  i.  e 
at  h = 1 and  T = 400  K. 

D = D^(T/T^)1-67/(p/pP 
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Schematic  of 
and 

TABLE  6.2 
Diffusion  Kinetic 
Guide  to  Figures 

Systems 

h,  T & T 

a 

slow 

border 

fast 

h=l 

T =400  K 

T,  =.004  s 
d 

1=1 

T =1.5  s 
c 

i=10 

T =.0015  s 
c 

6.2,6.3,6.18 

6.4,6.5,6.19 

h=3067 
T=400  K 
T =1.3(-6)  s 

i=.0367 

T =1.3(-5) 
c 

i=3.67 

s =1.3(-6)  s 

1=367. 

T =l.3(-7)  s 
c 

6.6,6.7,6.20 

6.8,6.9,6.21 

6.10,6.11,6.22 

i=1.44(-9)  i=1.44(-7)  i=1.44(-5) 

\ =l-2(-8)  s x^=l.2(-9)  s =1.2(-10)  s 


h = 3.34(6) 
T=300  K 
T =1.2(-9)  s 


6.12,6.13,6.23  6.14,6.15,6.24 


6.16,6.17,6.25 
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Figure  6.1.  Medium  density-source  strenth  grid. 

Points  at  which  diffusion  kinetic 
runs  were  made  are  indicated. 
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Comparison  of  Time  Courses . 

Comparing  first  of  ail  three  time  courses  for  slow  chemistry 
(t^  - IOt^)  we  note  that  in  general  the  majority  product  species — O2  or 
crosses  the  lower  range  of  concentration  levels  of  the 
initial  species  (i.e.  for  H2  and  0)  at  the  last  time  on  the  graph,  well 
after  the  concentrations  of  ail  these  initial  yield  species  have 
converged  as  a result  of  diffusion. 

Looking  next  at  the  borderline  cases  (1  atm  and  liquid  density),  we 
see  that  for  the  1 atm  case  OH  and  to  a lesser  extent  H shows  some 
effect  of  chemistry  in  their  time  courses,  while  both  show  about  equal 
effect  for  the  liquid  density  case. 

In  the  low  density,  somewhat  fast  chemistry  case,  it  can  be  seen 
that  OH  and  0 both  show  a strong  effect,  and  H2O2  and  HO2  also  by 
transmission  of  effect  from  OH  and  0 (OH  is  precursor  of  H2O2;  OH  and  0 
are  precursors  of  O2  which  reacts  with  H to  yield  HO2). 

We  may  next  compare  the  borderline  and  fast  chemistry  cases  with 
the  corresponding  pure  kinetic  cases  (Figs.  5.2a,  e and  i). 

1.  In  the  low  density  case  the  difference  diffusion  makes  is  that 
the  fast  diffusing  species,  H and  H2,  level  off. 

2.  In  the  1 atm  case  and  liquid  density  cases,  the  fast  chemistry 
case  was  chosen  to  be  faster  than  in  the  low  density  case 

^^c  " opposed  to  = t^/2.5).  The  choice  of  a 

lower  injection  rate  in  the  low  density  case  was  not  completely 
arbitrary.  Looking  to  the  medium  density-source  strength  grid 
(Fig.  6.1),  one  notes  that  this  injection  rate  is  already  near 
the  depletion  limit.  For  these  fast  chemistry  cases,  the  time 
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courses,  as  one  would  expect,  continue  for  a longer  time  to 
parailel  the  corresponding  pure  kinetic  time  courses. 


Comparison  of  Radial  Profiles. 

The  slow  chemistry  profiles  are  given  for  the  sake  of  comparison. 
Otherwise  they  provide  no  intrinsically  new  Information.  Let  us  first 
consider  the  fast  chemistry  case  at  low  density.  The  H profiies  are 
essentially  pure  diffusion  profiles.  H is  unreactive  at  this  low 
density,  as  noted.  The  profiies  of  OH  and  0 are  almost  identical  to 
each  other,  other  than  for  a scaling  factor.  This  reflects  the  role  of 
reaction  3.  Thus  0 present  at  later  times  becomes  primarily  more  and 
more  0 split  off  from  OH,  rather  than  initial  yield  0.  It  would  be 
interesting  to  calculate  age  distributions. 

We  next  consider  the  profiles  of  H2O2,  HO2  and  O2,  referring  to 
Figs.  5.13  e,  f and  g for  comparison.  Let  us  consider  first  the 
profiles  of  H2O2  and  O2.  The  reactant  species  which  produce  them  (OH;  0 
and  OH)  have  essentially  conformal  profiles.  They  themselves  have 
nearly  identical  diffusion  coefficient.  Their  profiles  differ  in  shape, 
that  of  O2  reflecting  an  inert  species,  which  it  is,  that  of  H2O2 
indicating  some  degree  of  chemical  destruction,  reflecting  the  fact  that 
H2O2  is  scavenged  by  OH  and  H (reactions  5 and  14  — see  Fig.  5.13g). 

We  may  compare  the  borderline  and  fast  chemistry  profiles  at  1 atm 
with  the  corresponding  profiles  at  1 atm  in  Chapter  III  for  the  isolated 
reaction,  H + OH  + M — > H2O  + M.  There  are  but  only  slight  differences 
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in  the  H and  OH  levels.  Referring  to  Fig.  5.14  a and  b,  we  note  that 

this  is  already  the  dominant  reaction  for  both  H and  OH  at  this  density 
already. 

We  next  consider  the  liquid  case.  Here  the  borderline  case  is  also 
not  extremely  interesting.  Going  on  to  consider  the  fast  chemistry 
case,  we  observe  that  the  depression  of  the  OH  peak  in  comparison  with 
the  H peak  is  even  greater  here  than  in  the  1 atm  case. 

In  contradistinction  to  the  low  density  and  1 atm  case,  in  the 
liquid  density  case  0 is  slow  to  react,  and  its  profiles  look 
essentially  like  pure  diffusion  contours.  Looking  at  Fig.  5.19d,  we 
note  only  one  destruction  reaction  (15),  which  offsets  only  one  third  of 
the  source  term  for  0.  As  one  can  see  from  Fig.  5.18  (reaction  rates), 
not  much  else  is  happening  in  the  liquid  density  case,  other  than  for  H 
+ OH  and  a brief  transient  in  OH  + OH  recombination.  That  the 
radj-olysis  of  liquid  water  by  an  electron  beam  is  very  slight  was 
observed  historically  by  Lea  (1946,  1947). 


Comparison  of  Near  Steady  State  Results 

It  was  noted  in  Chapter  III — in  regard  to  the  simple  diffusion 
kinetic  system  of  the  single  reaction,  H + OH,  that  the  most  physically 
realistic  or  depictive  presentation  of  the  results  is  given  by  a time 
sequence  of  radial  profiles  of  all  species  together.  Each  of  these 
plots  is  a snapshot  of  the  system  at  a specific  time.  For  the  case  of  a 
continuous,  steady  electron  beam,  this  time  sequence  would  depict  the 
approach  to  steady  state.  Thus,  in  this  case,  the  steady  state  snapshot 
is  of  particular  interest. 
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Because  it  becomes  expensive  to  run  the  code  out  to  extremely  long 
times,  it  is  necessary  to  settle  on  near-steady  state  results.  This  is 
a small  detraction.  For  only-diffusing  species,  the  time  course  is 
known,  anyway.  In  the  case  of  fast  chemistry  the  convergence  on  steady 
state  IS  abrupt  in  the  core  and  earlier  in  the  limb  than  in  the  only- 
diffusing  case.  The  treatment  in  Chapter  III  is  indicative. 

Snapshots  of  near  steady  state  at  the  last  time  calculated  for  in 
all  eight  cases  are  shown  in  Figs.  6.18-25.  In  the  very  slow  chemistry 
case  (Fig.  6.18),  the  full  steady  state  profiles  for  the  only-diffusing 
species  (the  'initial  yield'  species)  are  overlaid. 

In  order  to  facilitate  more  meaningful  intercomparison  between 
densities,  it  would  have  been  advantageous  to  have  implemented  equal 
reduced  final  times  (all  cases  at  a given  density  were  taken  out  to  the 
same  final  time).  Something  close  to  this  was  done.  Noting  that 
t'  = D3^d^  = T)yd\t  = 250ht  (using  = 4000  cm^/s),  the  reduced 
final  times  for  the  three  densities,  h = 1,  h = 3067,  and  h = 5.56x10^ 
(liquid),  are  250,  268.  and  150.  respectively.  Since  restart  diskflles 
were  recorded  to  tape  for  all  runs  made,  the  runs  at  liquid  density  can 
be  continued  to  t ' = 250,  should  it  be  desired. 

By  means  of  these  snapshots  of  near  steady  state  profiles,  it  is 
easier  to  make  qualitative  comparisons,  for  example,  regarding  level  of 
H2,  H2O2  or  O2  production,  comparison  of  those  levels  with  the  final 
ambient  levels  of  the  radical  species,  H,  OH,  0 and  HO2— both  with 
respect  to  densities  and  with  respect  to  injection  rate  at  a given 
density,  and  so  on. 
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At  a quick  glance,  we  note  the  following  significant  features. 

(1)  For  both  the  1 atm  cases  and  the  liquid  density  cases,  OH  can  be 
seen  to  fall  off  more  steeply  in  the  borderline  case  than  in  the  late 
chemistry  case.  In  the  early  chemistry  case  this  fall  off  becomes  still 
more  steep.  The  relative  peak  level  is  significantly  decreased  as  well. 

(2)  In  the  borderline  and  fast  chemistry  cases,  the  stable  or  relatively 

stable  products  and  the  radicals  are  generally  mixed  together  in  order 

of  magnitude  range.  No  simple  formula  differentiates  them.  In  the  fast 

chemistry  cases  (at  1 atm  and  for  liquid  density),  OH  and  HO2  are 

present  at  levels  significantly  lower  than  for  any  of  the  other 

species.  In  the  case  of  OH,  this  reflects  its  high  reactivity.  For 

HO2,  it  reflects  its  low  production  rate.  (3)  The  peak  steady  state 

concentration  levels  for  both  the  1 atm  fast  chemistry  case  and  the 

liquid  density  fast  chemistry  case  are  both  in  the  same  order  of 

magnitude  range  (10  ^ - 10  ^ mole/cc).  The  reason  for  this  is  as 

follows.  Both  are  'fast  chemistry'  cases,  i.  e.,  x = 0.1  x . 

c d 

Recalling  the  formulas  for  x and  x,,  this  relation  between  x and  x 

c a c d 

implies  the  following  relation  between  g and  h:  g/h  = 250.  Noting  that 

the  final  levels  are  roughly  proportional  to  gx^,  since  levelling  off 
occurs  right  around  or  not  too  long  after,  it  follows  that  these 
levels  are  proportional  to  (g/h)^-^^. 

Other  systeraatics  could  be  gleaned  from  these  snapshots  with  more 
effort.  However,  many  of  the  more  interesting  aspects  of  these 
diffusion  kinetic  systems  are  time  dependent  in  nature,  and  can  only  be 
apprehended  in  that  context.  These  include  (1)  the  pure  kinetic  basis 
for  comparison  with  the  full  diffusion  kinetic  case,  the  analysis  in 
the  previous  paragraph  in  terms  of  fundamental  time  dependences,  namely. 
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Figure  6.3.  Time  sequence  of  radial  profiles 
for  low  density  gas,  slow 
chemistry  case  for  (c)  H0„ . 
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Figure 


6,4.  Time  course  at  r = 0 for  T = 400  K 

E = 1 keV,  n~=  6 x lO^^  ///^c  and  I = 10^^  E/sec 
(low  density  borderline  chemistry)  . 
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Figure  6.6.  Time  course  at  r = 0 for  T = 400  K, 

E = 1 keV,  n = 1.84  x lO^^  ///cc  and 
I = 3.67  X 10^*^  #/sec  (1  atm  gas,  slow 
chemistry) . 
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Figure  6.7.  Time  sequence  of  radial  profiles 
for  1 atm  gas,  slow  chemistry 
case  for  (g)  HO.,. 
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Figure  6,8.  Time  course  at  r = 0 for  T = 400  K. 

= 1 keV,  n = 1.84  x lO^^  ///cc  and 
I = 3.67  X 10^  #/sec.  (1  atm  gas,  borderline 
chemistry) . 
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Fxgure  6.9.  Time  sequence  of  radial  profiles  for  1 atm  gas,  borderline  chemistry  case  for 
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Figure  6. 10. 


Time  course  at  r = 0 for  T = 400  K, 

~ 1 keV,  n = 1.84  x 10  if/cc  and 
I ~ 3.67  X 10^  iU sec  (1  atm  gas,  fast 
chemistry) . 
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Figure  6.11.  Time  sequence  of  radial  profiles  for  1 atm  gas,  fast  chemistry  case  for 
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Figure  6.11.  Time  sequence  of  radial  profiles  for  1 atm  gas,  fast  chemistry  gas  case  for 
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E = 1 keV,  n = 3.34  x 10^^  #/cc  and 
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Figure  6. 14.  Diffusion  kinetic  time  course  at  r = Q for 
T = 300  K,  = 1 keV,  n = 3.34  x 10^2 
and  I = 1.44  x IQ^  ///sec  (liquid  density, 
borderline  chemistry) . 


229 


o 


230 


(oo/afoui)  t uoT^Barjiraouoo  i 
o c 


2 ° 

■-I  iH 

(oo/3Xom)  uoT^Bj^^-uaouoo 


o 

CN 


o 

00 


e 

c 


o 

sr 


Figure  6. 15.  Time  sequence  of  radial  profiles  for  liquid  density,  borderline  chemistry  case  for 
(c)H„,  (d)  0. 
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Figure  6.16.  Time  course  at  r = 0 for  T = 300  K 

E,  = 1 keV,  n = 3.34  X lo22  L 
I - 1.44  X 10'  #/sec  (liquid  density, 
fast  chemistry). 
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Figure  6.17.  Time  sequence  of  radial  profiles  for  liquid  density,  fast  chemistry~case  for 
(a)  H,  (b)  OH. 
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Figure  6.17.  Time  sequence  of  radial  profiles  for  liquid  density,  fast  chemistry  case  for 
(e)  H0„,  (f)  0.. 
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Figure  6.18.  Snapshot  of  near  steady  state  profiles  for 
all  species  for  low  density,  slow  chemistry 
case,  t = 1 sec.  steady  state  pure  diffusion 
limit  overlaid  for  only-diffusing  species. 
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Figure  6.19.  Snapshot  of  near  steady  state  profiles  for 
all  species  for  low  density,  borderline 
chemistry  case,  t = 1 sec. 
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Figure  6.20.  Snapshot  of  near  steady  state  profiles 

for  all  species  for, 1 atm,  slow  chemistry 
case,  t = 3.5  X 10  sec. 
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Figure  6.21.  Snapshot  of  near  steady  state  profiles 
for  all  species  for  1 atm,  borderline 
chemistry  case,  t = 3.5  x 10~^  sec. 
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Figure  6.22.  Snapshot  of  near  steady  state  profiles  for 
all  species  for  1 atm,  fast  chemistry  case, 
t = 3.5  X 10  sec. 
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Snapshot  for  near  steady  state  profiles  for 
all  species  for  liouid  density,  slow  chemistry 
case,  t = 1.8  X 10  ' sec. 


Figure  6.23 
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Figure  6.24.  Snapshot  of  near  steady  state  profiles  for 
all  species  for  liquid  density,  borderline 
chemistry  case,  t = 1.8  x 10~' 


sec . 
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Snapshot  of  near  steady  state  profiles  for 
all  species  for  liquid  density,  fast  chemistry 
case,  t = 1.8  X 10"'  sec. 


Figure  6.25. 


246 


reaction  time  and  diffusion  time;  (2)  fundamental  diffusion  kinetic 
results  for  single  body  decay  and  for  both  one  and  two  species  two  body 
decay,  presented  in  Chapter  III;  (3)  time  dependent  diffusion  results, 
for  example,  diffusion  time,  but  more  extensively,  the  Green's  function 
approach  in  the  late  chemistry  case;  (4)  and  so  on. 

In  summary,  it  seems  fair  to  note  that  this  alternate  mode  of 
presentation  of  the  results  has  enabled  us  to  draw  one  very  significant 
conclusion  which  would  probably  have  escaped  us  otherwise.  This  mode  can 

properly  be  said  to  complement  the  modes  of  presentation  implemented  up 
to  this  point. 


Energy  Dependence 

It  would  be  useful  to  form  an  idea  of  the  dependence  of  the 
diffusion  kinetics  trends  in  g and  h,  equivalently  in  i and  h,  on 
..ncident  energy,  E^.  In  Chapter  II  we  concluded  that  the  spatial  yields 
could  be  represented  by  the  form 


y.(p,z)  = A,  (E  ) e 
i io,s  o 


-3. 


(6.1) 


where 


ij^^o^  ^ /100)I  (6.2) 


ij 


and 


5 = /Y^E^)p2  + (z  - z^(E^))2  - e(E^)(z  - z^(E^)) 


(6.3) 
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where  z (E  ) = z'(E  )/d. 
o o o o 


It  was  also  found  that  a reasonably  equivalent  concentric  form 
(i.e.,  one  for  which  e = 0)  could  be  constructed  by  making  the 
transformations 


e > e =0 
con 


Y > Y = Y 
con 


8 ^ 8^^^  = e(l  - e2)2/3 


(6.4) 


and  adjusting  by  nonlinear  least  squares  fit.  From  this  concentric 

form  an  equivalent  spherical  form  is  found  at  once  by  making  the 
transformations 


e 

con 


0 ->•  e , = 0 
sph 


Y ->■  Y = 1 

con  sph  . 

e ->■  s ^ = Y ^ e 

con  "^sph  ' ^con 


(6.5) 

z'  ^ Z*  ^ = Z* 

o,con  OjSph  o,con 


The  form  of  y(p,z)  now  simplifies  to 


y^(r)  = 


A.  (E  )e 
io,s  o 


-6, 


(6.6) 


where  r = /p  + (z  - z ). 

o 


248 


The  question  of  what  to  do  with  this  double  exponential  form, 
considering  that  the  development  of  all  of  our  diffusion  kinetic  results 
so  far  has  been  based  on  a single  exponential  form  for  the  source  terms, 
stands  now  in  our  way.  In  either  extreme,  where  one  of  the  two 
exponentials,  the  short  range  or  the  long  range,  vastly  predominates 
over  the  other,  the  problem  is  solved — the  source  term  effectively 
degenerates  to  a single  exponential.  The  ratio,  J approaches 
a maximum  value  of  4.4  (recall  k = 4.4  in  Eq.  (2.29))  as  decreases, 
and  approaches  a minimum  value  of  1 as  E^  increases.  The  ratio 
^io,s^®i  approaches  a maximum  value  of  1 as  E^  decreases,  and  a 
minimum  value  of  about  0.1  as  E increases.  For  E > IkeV  the  short 
range  term  can  be  neglected,  while  for  E^  < 300  eV  the  long  range  term 
can  be  neglected.  In  between,  we  may  approximate  matters  by  a single 
exponential  for  which  B is  assigned  some  intermediate  value  between  the 
two  extreme  limits.  Since  this  assignment  is  an  estimate,  a confidence 
interval  should  be  assigned  along  with  it. 

Ibe  interesting  quantities  whose  dependence  on  E 

o we  would  like  to 

identify  are  and  t^.  Let  us  first  consider  t^.  The  diffusion 

time,  T^,  is  given  by  by  l/D^d^B^(E^)h.  It  is  proportional  to  the 

square  of  the  range  of  the  deposition,  and  therefore  increases  with 

increasing  E^.  In  the  Bom  limit,  it  therefore  goes  as  E The 

o 

reaction  time,  r^,  is  given  by 


T 

C 


- go(E^))/100)E^(B^E^)d^I/8T;)  , d < 4.6(-4)  gm/cc 


d > 4.6(-4)  gm/cc 


^/^‘^6b^SoH(E^)/100)E^(B^(E^)d3l/8ir)Cj^  ^ 


(6.7) 
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or,  equivalently  in  terms  of  g and  h,  by 


/ 


T 


C 


For  > 200  eV  the  g-values  are  nearly  constant.  dependence  comes 

in  through  the  factor  of  E and  from  B(E  ).  In  the  Born  limit,  it 

u o ’ 


increase  much  more  slowly  with  E^  than  the  diffusion  time.  All  other 
things  being  equal,  the  chemistry  becomes  progressively  earlier,  in  fact 
much  earlier,  as  the  incident  energy  of  the  electron  beam  is  increased. 

For  very  high  incident  energy,  the  contours  of  deposition  become 
more  and  more  singly  confocal,  so  that  a concentric  approximation,  let 
alone  a spherical  one,  becomes  less  and  less  adequate.  The  problem  then 
becomes  irreducibly  two  dimensional. 

In  the  extreme  limit,  the  geometry  of  the  problem  approaches  that 
of  a cylinder,  which  is  also  one-dimensional — specifically,  a a cylinder 
with  a bulb  on  the  end  where  the  electron  and  its  high  energy 
secondaries  have  significantly  slowed  down,  and  are  giving  rise  to  many 
low  energy  secondaries.  At  not  such  extreme  high  energy  limit 
(relativistic),  this  bulb  engulfs  the  narrow  core  or  stem  of  higher 
energy  electrons  which  give  rise  to  it,  resulting  in  the  contours  we 
have  calculated  at  10  keV  (see  Fig.  2.10),  and  which  Griin  (1956)  and 


therefore  goes  as  E 


2.1 


The  reaction  time  is  therefore  seen  to 


o 
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Cohn  and  Caledonia  (1970)  report  at  incident  energies  of  from  5-54  keV 
and  2-5  keV,  respectively.  In  Chapter  II  we  have  referred  to  this 
configuration  as  'penumbra’  and  'core.' 

Let  us  consider  for  a moment  such  a plume  shape  in  detail  before 
gOj.ng  on.  It  is  characterized  by  a forward  and  a backward  fail  off, 
e^CE^)  = e(E^)(l  - e(E^))  and  B^(E^)  = 6(E^)(1  + e(E^)).  Because 
~ — £ (E^)  for  high  E^,  the  fitted  contours  are  nearly 

Cj.rcular  (this  is  a simplified  approximation  in  the  central  core)  about 
a common  'affine  focus'.  The  geometry  of  the  backplume  is  something  in 
between  a nearly  flat  lens-like  shape  and  a hemisphere  half-centered  on 
the  affj.ne  focus.'  The  forward  plume  resembles  a somewhat  deformed 
sphere. 

The  spherlcalization  of  the  entire  form  (see  above)  more  strongly 
resembles  the  forward  plume.  This  reflects  the  relative  sizes  of  these 
two  regions.  From  Eqs.  (6.4)  and  (6.5)  it  follows  that  the 
sphericalized  fail  off,  is  equal  to  B(1  - e^),  which  lies  in 

between  the  forward  fall  off,  6(1  - e),  and  the  backward  fall  off, 

6(1  + e),  but  closer  to  the  forward. 

This  essay  in  estimation  has  been  presented  in  order  to  lay  a 
foundation  for  a future  link-up  of  the  spherical  form  source  term 

^iii®tic  results  obtained  here  and  future  work  with  source 
terms  of  more  general  shape.  Such  future  work  is  seen  as  a possibility 
via  an  upgrading  of  DIFKIN  to  calculate  for  a generalized  polar  as  well 
as  generalized  radial  grid.  Specifically,  the  affine  singly  confocal 
coordinate  system,  or  even  a more  general  radial-type  coordinate  system, 
is  seen  as  the  most  useful  route  to  take  for  this  purpose. 
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Heating  Effects 

In  the  radiological  problem  addressed  in  this  work  certain  physical 
assumptions  have  been  made.  The  first  of  these  is  that  the  system 
remains  dilute,  in  other  words,  that  a steady,  uniform  background  of  H2O 
be  maintained.  This  is  necessary,  in  the  first  place,  because  the 
calculation  of  the  spatial  yields  is  based  upon  this  assumption.  It  is 
necessary,  in  the  second  place,  because,  should  it  not  be  the  case, 
then,  because  of  the  law  of  partial  pressures,  a convective  regime  will 
be  set  up,  which  contradicts  the  assumption  of  the  diffusion  kinetic 
equations,  that  the  medium  is  stagnant.  The  "depletion  time,"  defined 
in  Chapter  V,  can  serve  to  identify  an  upper  limit  on  injection  rate 
mandated  by  this  requirement.  Of  the  eight  hypothetical  cases 
considered,  one,  the  high  injection  rate  case  for  low  density  gas  lies 
dangerously  near  this  limit.  In  this  case,  H2O  level  is  down  to  94.6% 
r— 0 at  0.1  sec,  and  92.8%  at  1 sec.  In  all  the  other  cases,  H2O 
depletion  is  negligible. 

A yet  more  severe  limit,  the  upper  limit  on  injection  rate  due  to 
heating,  must  next  be  addressed.  A full  treatment  of  this  problem  would 
have  to  take  into  account  the  energetics  of  the  chemical  reactions.  In 
lieu  of  this,  a rough  order  of  magnitude  estimate  is  made  here,  arguing 
from  conservation  of  energy  and  the  equipartition  theorem. 

The  time  rate  of  energy  delivered  to  the  system,  or  power  input  to 
the  system,  is  given  by  P = lE^.  In  lieu  of  a more  involved  statistical 
mechanical  derivation,  as  an  order  of  magnitude  estimate  we  may  set 

p(r)  = IE^(8^d^/8Tr)e 


(6.9) 
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Divj-ding  by  7nk/2,  where  n,  as  used  before,  is  the  number  density 
(assigning  3kT/2  for  kinetic  energy,  kT  for  vibrational,  and  kT  for 
rotational  energy),  one  obtains  a spatially  distributed  heat  source, 
y^-^Cr),  as  follows: 


yj.j^(r)  = IE^(S^d^/87r)e"^'^''/(7nk/2) 


T a 2 

^o°o  0 

■ ('28.N  k " 
o 


0.4  CK/sec)  (E^/1  keV)[B(E  )/B(l  keV)]^.h^i 


(6.10) 


Here  ^ is  the  gram-molecular  weight  of  H2O. 

Taking  into  account  thermal  conductivity  of  the  medium,  the  heat 
equation,  or  thermal  diffusion  equation. 


at 


(6.11) 


follows,  where  c is  molar  specific  heat. 


Cy  Q is  molar  density. 


temperature,  and  temperature  dependence  of  the  thermal 
conductivity,  K,  has  been  neglected.  Noting  that  K = K^(T/T^)*^  where 
q - 1.83  and  = 273  K and  = 3.61x10  ^ cal/gm/cm/sec/°K  for  steam 
(Kanury,  1975),  and  recalling  the  formula  for  diffusion  time  from 
Chapter  III,  we  obtain 
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^th  ~ (T/T  )‘^/h  sec  (vapor)  (6.12) 

H^O  2 

The  equivalent  diffusion  coefficient  differs  on  the  low  side  of  the 
diffusion  coefficient  of  H2O  due  to  effects  of  nonmonatomic  gas.  For 
liquid  water,  the  still  lower  value,  K = .00146  cal/gra/cin/sec/°K  at 
300  K (CRC  Handbook,  1975),  reflecting  collective  effects. 

Recalling  the  definition  of  diffusion  time,  it  follows  that 
= T(r=0,«>)  is  given  by 


Tf  = yth^^=o>^th 

= (0.4h^i)(.0019/h)  °K 


= 0.00076hi  (°K). 


(6.13) 


The  values  of  T^  which  follow  from  Eq.  (6.13)  for  the  eight  cases 
considered  are  given  in  Table  6.3.  These  values  should  only  be  taken  as 
order  of  magnitude  estimates,  and  could  be  wrong  on  either  the  low  or 
the  high  side.  All  the  same,  we  are  led  to  conclude  that  two  and 
possibly  three  of  the  cases  considered  are  physically  ruled  out  as  a 
result  of  overheating.  (In  these  two  cases,  the  significance  of  the 
estimate  is  only  that  the  temperature  Increase  is  too  large.  For  such  a 
strong  heating  effect,  the  T dependence  of  K,  along  with  convective 
effects  of  thermal  expansion,  invalidate  Eq.  (6.11).)  Still,  something 
has  been  gained  by  considering  them  insofar  as  our  understanding  of  the 
behavior  of  diffusion  kinetic  systems  in  general.  For  example,  the 
character  of  results  for  the  isolated  H + OH 


system  considered  in 
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TABLE  6.3 

Rough  Estimate  of  Heating  Effects  of  Electron  Beam 


case 

n I 

(#/cc)  (///see) 

h 

i 

g est 

Tf(r=0) 

(°K) 

1 

6.(15) 

1.(12) 

1 

1 

1 

0.0008 

2 

6.(15) 

1.(18) 

1 

1.(6) 

1.(6) 

800. 

3 

1.84(19) 

3.67(10) 

3067. 

0.0367 

1.06(9) 

0.045 

4 

1.84(19) 

3.67(12) 

3067. 

3.67 

1.06(11) 

4.5 

5 

1.84(19) 

3.67(14) 

3067. 

367. 

1.06(13) 

450. 

6 

3.34(22) 

1.44(3) 

5.56(6) 

1.44(-9) 

2.48(11) 

4. (-7) 

7 

3.34(22) 

1.44(5) 

5.56(6) 

1.44(-7) 

2.48(13) 

4. (-5) 

8 

3.34(22) 

1.44(7) 

5.56(6) 

1.44(-5) 

2.48(15) 

0.004 

h = d/do  = where  n = 6.(15)  #/cc 

i = I/Iq  where  = 1.(12)  #/sec 

g = ih-^ 

Tf(r=0)  = Yth(’^=0)T^^  = (.4g/h)(.002/h)  = .0008g/h^ 

= (.4ih^)(.002/h)  = .OOOSih 

for  ambient  T = 400  °K,  vapor 

= (*4g/h)(.000094/h)  = .000038g/h 
= (.4ih^)(.000094/h)  = .000038ih 

= .122(-17)g 
= 209i 

for  ambient  T = 300  °K,  liquid  (h  = 5.56(6)) 
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Chapter  III  remain  essentially  unchanged  in  going  over  to  the  physically 
more  realistic  system  in  the  liquid  state — only  spatial  and  temporal 
scales  would  be  changed. 


CHAPTER  VII 

DISCUSSION  AND  CONCLUSION 

Comparison  with  Results  of  the  Steady  State  Code,  DIFKT 
The  work  embodied  in  this  thesis  was  originally  envisaged  as  a 
backup  approach  to  a parallel  effort  to  calculate  steady  state  diffusion 
kinetics  for  the  same  inputs,  namely,  the  Ph.  D.  thesis  work  of  Rio 
(1983).  That  work  however  remains  in  a comparatively  rough  state  of 
development.  As  a result  of  continuation  of  this  steady  state  diffusion 
kinetics  work  by  collaborative  effort  within  this  group.  Dr. 

Dayashankar,  Professor  Green  and  the  author,  certain  numerical 
difficulties  of  the  code  have  come  to  light.  These  problems  include; 

4 

(1)  In  assorted  checks  devised  to  test  the  accuracy  of  this  code 
for  simple  problems  of  pure  diffusion,  diffusion  kinetics  of 
^ic^ic  body  decay  of  one  species,  and  of  two  body  decay  of  one 
or  two  species,  the  code  has  shown  unacceptable  levels  of 
inaccuracy. 

(2)  The  code,  based  on  the  finite  element  partial  differential 
equation  solving  package,  TWODEPEP  (IMSL,  1980)  requires  an 
i^-^tisl  distribution  of  nodes,  or  vertices,  of  the  finite 
elements.  If  these  nodes  are  not  aligned  with  the  contours  of 
the  source  terms,  they  should  be  expected  to  migrate  to  where 
they  are  needed.  Apparently,  this  is  not  happening,  or  if  it 
is,  something  else  is  wrong  in  its  functioning  as  evidenced  by 
essentially  random  patterns  in  the  alignment  of  the  peaks  of 
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the  calculated  concentration  contours  with  the  peaks  of  the 
source  terms.  Specifically,  the  peaks  in  the  concentrations  are 
showing  up  more  or  less  around  the  origin  (where  the  peak 
initial  node  density  is  located)  or  between  it  and  the  axial 
location  of  the  source  term  peaks,  but  still  hanging  close  to 
the  origin. 

(3)  A set  of  strange  systematics  obtains  when  ratios  are  taken 
between  the  peak  values  obtained  by  this  code  and  those 
obtained  by  DIFKIN.  These  systematics  may  hold  the  key  to  the 
-i.nner  working  of  TWODEPEP  in  relation  to  its  implementation  in 
DIFKT,  and  the  problems  which  are  being  encountered  in  the 
functioning  of  DIFKT. 

(4)  The  diffusion  kinetic  calculation  being  attempted  with  DIFKT 
for  the  deposition  problem  is  for  the  full  kinetic  system  of 
i^sutrals,  ions  and  thermal  electron.  The  species  which  show 
themselves  to  be  most  sensitive  to  the  nonalignment  of  node 
distribution  and  source  contours  are  the  ionic  species  (see 
Figs.  7.1-2).  Their  concentrations  show  a strange  but 
consistent  asymmetry  as  a result.  The  vastly  different  time 
scale  of  these  reactions  from  that  of  the  neutrals  (for  the 
conditions  being  calculated)  presumably  plays  a role.  The 
global  recombination  approach  devised  here  may  alleviate  this 
problem  for  DIFKT. 

In  results  of  runs  of  DIFKT  at  100  eV,  where  the  spatial  grid  was 
left  the  same  as  for  1 keV,  as  a result,  the  discretization  of  the 
source  was  crudified  by  about  a factor  of  10  (since  at  this  low  energy 
the  short  range  term  predominates  over  the  long  range,  so  the  decrease 
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Figure  7.2.  Results  of  DIFKT  for  the  low  density,  borderline  chemistry  case,  (unpublished 
work  of  Rio,  Dayashankar,  Schippnick,  and  Green)  for  (a)  neutral  species, 
and  (b)  ionic  species  and  thermal  electron. 
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in  range  is  still  more  than  what  it  would  be  if  the  dependence  of  the 
long  range  were  the  only  factor).  As  a result,  in  these  runs  the  effect 
of  the  axis  offset  discrepancy  with  the  location  of  the  node 
distribution  peak  (the  origin)  was  nil.  In  this  case,  the  concentration 
contours  of  the  ions  did  not  show  an  asymmetry.  In  fact,  they  presented 
extremely  narrower  profiles  in  this  case.  This  result  is  an  interesting 
one,  and  may  provide  a further  clue  in  our  search  for  what  is  going 
wrong  in  the  functioning  of  this  code. 

One  aspect  of  the  functioning  of  DIFKT  which  may  play  a role  in  the 
above  problems  is  the  mode  of  coupling  it  implements.  As  a result  of  the 
limitation  of  TWODEPEP  to  couple  five  separate  variables  at  one  time, 
the  coupling  between  all  concentrations  is  realized  by  coupling  them  in 
fixed  groups  of  five  per  iteration;  the  coupling  between  groups  is 
effected  by  the  iteration  step  in  the  reassignment  of  values.  The 
question  naturally  arises— what  effect  is  there  in  grouping  one  way 
rather  than  another? 

DIFKIN  may  be  of  use  in  the  upgrading  of  DIFKT  in  the  following 
way.  DIFKT  allows  for  initial  estimates  of  the  solution.  Convenient 
analytic  forms  have  already  been  established  here  for  most  of  the  forms 
which  the  steady  state  concentration  profiles  take.  In  this  way,  say,  in 
application  to  particularly  nonspherical  spatial  yield  contours,  the 
spherical  approximation  solution  can  be  obtained  via  DIFKIN,  modified  to 
a nonspherical  shape,  and  used  as  an  initial  estimate  in  DIFKT.  It 
would  be  interesting  to  see  what  effect  this  would  have  on  the  rate  of 
convergence  of  iteration,  and  on  the  level  of  accuracy  of  the  final 
answer,  in  DIFKT. 
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In  spite  of  the  defects  in  the  results  obtained  from  DIFKT  so  far 
they  are  in  fairly  reasonable  agreement  with  the  resuits  of  DIFKIN. 
Specj.f icaliy , DIFKT  has  been  run  for  the  two  low  density  cases  presented 
in  Chapter  VI,  very  late  chemistry  and  somewhat  fast  chemistry.  These 
results  (Rio,  Dayashankar,  Schippnick  and  Green,  unpublished  work)  are 
displayed  in  Figs.  7.1-2.  The  results  for  the  late  chemistry  case  are 
in  better  than  order  of  magnitude  agreement  with  the  results  obtained 
here  (Fig.  6.18).  This  is  interesting  in  regard  to  the  chemical  product 
species,  in  regard  to  the  distinction  between  near  steady  state  and 
actual  steady  state.  It  is  quite  possible  that  the  iterative  procedure 
acts  like  a pseudotime  (in  a very  distorted  sense,  but  comparable). 

This  idea  was  proposed  by  Rio  (private  communication,  1983). 

The  results  for  the  fast  chemistry  case  are  also  in  better  than 
order  of  magnitude  agreement  with  the  results  obtained  here  (Fig. 

6.19).  Two  points  regarding  the  comparison  here  and  in  the  previous 
case  should  be  noted.  First,  the  spatial  yields  inputs  implemented  in 
these  steady  state  calculations  are  essentially  spherical  (axis  ratio, 
.97).  Thus,  the  two  results  can  be  compared  straight  across  in  this 
regard.  (The  worrisome  axis  offset  here  is  8.6  cm.)  Second,  DIFKT  is 
calculating  for  the  full  ion-neutral-electron  system.  The  results 
obtaining,  if  they  be  true  results,  should  for  that  reason  not  be  in 
complete  agreement  with  the  results  obtaining  from  DIFKIN.  The  degree 
of  this  difference  resulting  from  the  global  instantaneity  approximation 
remains  to  be  determined.  In  this  regard  it  should  be  noted  that  the 
spatial  yield  inputs  used  in  this  calculation  by  DIFKT  were  not 
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corrected  for  any  discrepancy  in  charge  balance.  This  should  be  looked 
after  in  future  work,  as  to  the  sensitivity  of  the  calculation  to  this 
factor. 

reassuring  that  the  initial  ion  species  appear  at  the  lowest 
concentrations,  reflecting  the  efficiency  of  reaction  with  H2O.  It  is 

interesting  that  in  both  cases  H'''(HoO)o  j . + 

' 2 ^3  predominates  over  H^O  and 

H"'’(H20)2.  It  is  not  clear  how  this  relates  to  the  fact  that  its 

reaction  rate  with  the  electron  is  the  fastest  of  the  three.  It  may 

relate  to  the  fact  that  it  is  the  last  one  to  be  produced,  although  in 

the  pure  kinetic  case  that  fact  correlates  with  it  being  the  one  in 

lowest  ambient  concentration.  This  is  an  open  question  since,  for  one 

thing,  the  overall  accuracy  of  this  calculation  is  uncertain. 

We  conclude  that  the  code,  DIFKT,  has  good  promise,  but  still 
suffers  from  some  nagging  problems,  which  have  yet  to  be  resolved.  With 
more  work,  these  problems  should  be  resolved.  It  is  worth  noting  in 
this  regard  that  the  IMSL  package,  TWODEPEP,  has  been  thoroughly  tested 
and  used  by  many  users  over  the  last  10  years  or  so,  successfully  and 
accurately  (for  example,  see  Mooney,  1984).  Because  DIFKT  still  gave 
assorted  errors  when  tested  with  simpler  systems,  such  as  those 
considered  in  Chapter  III,  we  are  led  to  conclude  that  there  are  errors 
in  the  way  in  which  it  implements  TWODEPEP. 


Proposed  Experimental  Tests 

Experimental  identification  of  'initial'  spatial  yields,  i.  e., 
spatial  yields  of  various  excited  states  created  by  the  energetic 
degradation  of  the  primary  electron  and  its  secondaries,  etc.  is 
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facilitated  by  the  emission  of  characteristic  wavelengths  associated 
with  these  excited  states  when  they  decay  to  lower  states  and,  finally, 
ground  state.  On  the  basis  of  line  of  sight  detection  of  these 
emissions,  the  spatial  distributions  of  the  excited  states  which  emit 
them  may  be  inferred  via  Abel  inversion  (Cohn  and  Caledonia,  1970). 

The  earliest  such  work,  and,  to  this  day,  definitive  in  its 
results,  was  the  work  of  Grun  (1956)  for  5-54  keV  electrons  incident  in 
air  at  700  Torr.  The  later  work  of  Cohn  and  Caledonia  (1970)  was  for  2- 
5 keV  electrons  Incident  in  N2  at  0.69-2.8  Torr.  In  both  cases  spatial 
contours  are  shown  very  similar  to  the  spatial  contours  of  Monte  Carlo 
generated  spatial  yield  spectra  presented  in  Fig.  2.10.  Later  work  by 
Barrett  and  Hays  (1975)  obtained  similar  results. 

In  order  to  identify  spatial  distributions  of  chemical  species  (i. 
e.,  ground  state  species)  in  the  diffusion  kinetic  evolution  of  the 
irradiated  medium,  a different  approach  is  necessary.  If  the  species  of 
interest  undergoes  a transition  lying  in  a detectable  wavelength  band, 
with  linear  absorption  characteristics  (i.  e.,  obeying  Beer’s  law)  or  at 
least  known  absorption  characteristics,  then  its  presence  can  be 
detected  via  absorption  spectroscopy.  Again  line  of  sight  detection  can 

be  utilized,  and  a variation  of  Abel  inversion  devised  to  infer  spatial 
distribution. 

In  the  case  of  electron  bombardment  of  H2O,  the  only  species  whose 
characteristic  spectra  lie  in  the  optical  or  ultraviolet  are  OH  (309.5 
nm),  H (Lyman  a , 121.6  nm)  and  0 and  possibly  HO2. 

In  this  connection,  the  pulsed  radiolysis  work  of  Hanrahan  and 
coworkers  (this  campus,  see  for  example  Ramirez,  Bera  and  Hanrahan, 

1984)  IS  of  interest.  A brief  description  of  their  work  follows.  A 
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Febetron  electron  beam  generator  (Bera  and  Hanrahan,  1984)  is  used  which 
typically  produces  a 0.6  MeV  beam  of  about  3 nanoseconds  duration.  The 
beam  energy  can  be  taken  down  to  about  0.3  MeV  but  not  much  lower 
because  of  the  parallel  charge-series  discharge  arrangement  of  the 
capacitor  bank  used  to  generate  the  beam.  The  medium  consists  typically 
of  about  600  Torr  buffer  Ar  gas,  variable  amounts  of  water  vapor  (in  the 
range  0.5-15  Torr),  and  reactive  substrate.  Concentration  of  OH  is 
monitored  over  a 10  psec  interval.  This  measurement  is  not  spatially 
resolved.  Presumably,  it  is  possible  to  do  so. 

Whether  the  possibility  exists  of  carrying  out  continuous  beam 
experiments  which  would  correspond  exactly  to  the  configurations 
calculated  for  in  this  thesis  is  uncertain.  It  is  a generally  held 
opinion  (Hanrahan,  private  communication,  1984)  that  continuous  beam 
generators  generally  do  not  exist.  Another  consideration  is  the 
interrelationship  between  beam  energy,  medium  density,  spatial  range, 
duration  of  pulse  and  desired  duration  of  pulse.  As  density  increases, 
the  time  of  near  steady  state  decreases  (diffusion  time).  Spatial  range 
also  decreases.  A minimum  range  is  defined  by  maximum  spatial 
resolution  of  detector.  Most  probably,  a fast  chemistry  system  which 
does  not  overheat  (at  least  several  atmospheres)  would  require  too  high 
a spatial  resolution  to  be  practicable.  The  load  of  the  increased 
pressure  on  differential  pumping  of  the  electron  gun  (for  example  see 
Griin,  1956)  also  has  limit  and  imposes  one.  That  which  it  imposes  is 
the  convective  disturbance  near  the  gun.  Thus  a spatial  range  too  small 
would  be  entirely  engulfed  in  this.  The  condensation  point  limit  of  H2O 
vapor  pressure  must  also  be  taken  into  account. 
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The  above  discussion  is  meant  only  to  be  very  preliminary. 
Experimental  realization  of  continuous  beam  systems  may  be  possible. 
However,  the  indications  are  that  it  is  not.  As  an  alternate  route  by 
which  the  radiological  work  carried  out  in  this  group  may  link  up  with 
experiment,  the  pulsed  radiolysis  work  of  Hanrahan  and  coworkers  offers 
itself  as  a real  possibility.  Towards  this  end  the  following  plan  of 
work  appears  feasible:  1.  Extend  the  spatial  yield  spectrum  to  higher 
(relativistic)  energies.  2.  Calculate  the  spatial  yield  spectrum  for 
mixtures  of  gases,  according  to  the  prescription  given  in  Singhal  and 
Green  (1982),  which  is  based  on  the  near  universality  of  the  spatial 
yield  spectra  for  all  gases.  In  this  way,  electron  yield,  Ar  yield,  and 
the  initial  yields  of  H2O  radiolysis  are  easily  arrived  at.  3.  Modify 
DIFKIN  to  calculate  in  cylindrical  geometry.  4.  Run  DIFKIN  for  an 
initial  pulse  rather  than  a continuous  pulse  beam.  5.  Calculate  pure 
diffusion  results  via  Green's  function  integral  for  cylindrical  geometry 
for  initial  pulse.  6.  Calculate  pure  kinetic  results  for  initial  pulse, 
also  over  g-h  grid,  to  enable  isolation  of  diffusion  effect  in 
comparison  of  diffusion  kinetic  result  with  pure  kinetic  result. 


Summary  and  Conclusion 

The  problem  of  diffusion  kinetics  for  continuous,  steady  electron 
beam  injection  in  water  vapor  has  been  addressed  in  a time  dependent 
context.  Certain  approximations  have  been  implemented  in  order  to 
^^cilitate  calculation.  Quantitative  limits  of  the  ranges  of 
applicability  of  these  approximations,  and  quantitative  estimates  of  the 
accuracy  of  these  approximations  have  either  been  made,  or  the  means 
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with  which  to  do  so  have  been  produced.  For  example,  in  the  case  of  the 
approximation  of  global  instantaneity  of  charge  recombination,  it  was 
observed  that  the  effect  of  the  approximation  is  limited  to  the  H 
radical  yield,  as  far  as  the  neutral  chemistry  is  concerned. 

The  other  significant  approximation  which  has  been  made  is  the 
sphericalization  of  the  deposition  contours.  Referring  to  the  steady 
state  axisymmetric  diffusion  result  calculated  in  Appendix  C,  and  to  the 
degree  of  eccentricity  of  the  deposition  contours,  it  is  clear  that 

been  lost  by  this  for  contours  which  are  not  highly  confocal. 

It  has  been  striven  here  to  treat  every  aspect  of  the  problem  in  as 
quantitative,  rigorous  and  thorough  a manner  as  possible.  Where  there 
have  been  uncertainties,  such  as  in  regards  to  the  stiffness  problem,  or 
possible  numerical  artifacts  of  the  radial  grid,  they  have  been  pointed 
out. 

Data  display.  Close  attention  has  been  paid  to  the  display  of 
data.  For  the  sake  of  clarity  and  convenient  visual  display,  several 
results  have  been  presented  in  more  than  one  format.  The  case  of  net 
production  rates  is  particularly  relevant  because  it  is  so  difficult  to 
estimate  a slope  from  data,  and  because  this  result  is  at  the  heart  of 
the  stiffness  problem.  This  result  illustrates  the  importance  of 
checking  results. 

In  the  preliminary  work  to  this  thesis,  presented  in  Chapter  II, 
display  of  data  has  also  played  an  important  role,  and  can  be  expected 
to  play  an  even  more  important  role  in  the  continuation  of  this  work  in 
pursuit  of  an  optimal  analytic  form  by  which  maximal  accuracy  can  be 
achieved  in  the  analytic  representation  of  the  spatial  yield  spectra. 

What  is  envisioned,  along  with  residual  plots,  are  discrete  bin  band- 
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contour  plots,  where  intervals  are  contoured,  by  which  we  can  actually 
see  the  data.  Without  visual  aids  of  this  sort,  identification  of 
patterns  and  trends  is  severely  hindered  by  the  tediousness  of  manual 
plotting  of  a few  representative  cases. 

Expense  of  the  code.  A point  which  has  played  a role  of  some 
importance  in  this  work  has  been  the  expense  of  the  code.  Because  the 
implicit  linearized  calculation  of  the  chemistry  in  the  BOAT  code  solves 
a system  of  n simultaneous  linear  equations  at  each  grid  point  (where  n 
is  the  number  of  species),  the  cost  of  the  code  is  primarily  a function 
of  the  number  of  species. 

The  size  of  the  chemistry  friendly"  time  step  also  bears  on  the 
cost  of  the  code,  but  less,  because  the  major  source  of  cost  is  in  the 
logarithmic  growth  of  time  step  at  long  times  which  is  a result  of  the 
logarithmic  evolution  of  the  grid  point  distribution.  This  evolution  is 
best  kept  logarithmic  to  insure  that  the  grid  point  redistribution 
adequately  cover  the  evolving  steady  state  core.  However,  it  can 
probably  be  made  to  move  along  faster  without  significant  loss  of 
accuracy.  The  advantage  here  can  probably  be  enhanced  by  replacing 
Simpson  rule  interpolation  by  cubic  interpolation.  This  would  be 
beneficial  for  accuracy  in  general.  It  should  also  be  possible  to 
freeze  out  the  core  as  it  evolves.  In  this  way  a significant 
economization  could  be  attained. 

However,  if  the  full  system  of  ions,  neutrals  and  thermal  electron 
were  to  be  calculated  for,  the  expense  due  to  "chemistry  friendly"  time 
step,  as  well  as  the  expense  due  to  the  increase  in  the  size  of  the 
system  of  simultaneous  linear  equations  being  solved  for,  would  add  up 
to  a significant  increase  in  expense. 
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Calculating  on  an  IBM/370  type  machine,  one  marching  step  requires 
0.08  sec  of  CPU  time  for  a 21  point  grid  for  the  full  neutral  system. 

The  number  of  steps  in  the  diffusion  kinetic  runs  presented  in  Chapter 
III  and  Chapter  VI  gives  an  idea  of  the  typical  expense  of  this  code. 

At  the  extreme  of  slow  chemistry,  the  low  density,  1 atm,  and  liquid 
runs  took  19320,  17523,  and  16574  steps,  respectively.  At  the  extreme  of 
fast  chemistry,  the  1 atm  and  liquid  runs  took  30645  and  28788  steps, 
respectively. 

Checking  of  code.  The  code  has  been  checked  out  in  detail,  and 
found  to  give  very  accurate  results  when  compared  with  known  analytic 
test  cases  (Chapter  III).  These  test  cases  include  pure  diffusion, 
checked  against  the  Green's  function  calculation,  pure  kinetics  for  one 
and  two  species  systems,  steady  state  diffusion  kinetics  for  single  body 
decay,  checked  against  analytic  results,  and  steady  state  diffusion 
kinetics  for  one  species  two  body  decay,  checked  against  analytic 
results  in  its  asymptotic  limits. 

Approximations.  Calculating  for  an  incident  energy  of  1 keV,  which 
was  taken  as  a convenient  representative  value  from  which  meaningful 
extrapolation  both  above  and  below  is  possible,  provided  a further 
convenience  in  that  for  the  spatial  yields  at  this  incident  energy,  the 
short  term  component  is  negligible,  and  a single  term  exponential 
representation  is  possible.  As  a result,  the  analytic  work  pertaining  to 
diffusion  and  diffusion  kinetics,  in  which  it  was  desired  to  use  a 

source  term  of  the  same  form  as  in  the  physical  problem,  was  greatly 
simplified. 

However,  whether  single  or  double,  the  exponential  form  is 
deficient  in  the  following  respect,  namely  that,  the  slope  is  not 
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continuous  at  the  peak.  It  is  easy  to  construct  such  a form.  Two 
possibilities  are  (1)  in  general,  the  mapping,  e~®^  -*■  where 

f (0)  - 0 and  f(r)  ~ r as  r > and  (2)  the  weighted  difference  of 
two  exponentials,  - 3e~“'')/( ct-6) . The  latter  enjoys  the 

advantage  of  a linear  sum  in  the  pure  diffusion  result  of  its  Individual 
terms.  However,  once  the  nonlinear  kinetic  term  enters  in,  it  possesses 
no  more  of  an  advantage  than  the  other. 

Why  this  point  can  matter  is  that  the  'reaction  time'  of  the 
diffusion  kinetic  system  has  been  defined  in  terms  of  the  peak  values  of 
the  spatial  yields.  We  would  therefore  like  to  have  the  most  accurate 
representation  of  the  spatial  yields  that  we  can  obtain.  However, 
considering  that  the  overestimation  which  comes  as  a result  seems  to  be 
small,  considering  the  profiles  of  spatial  yield  spectra  and  deposition 
spectra  displayed  in  Chapter  II,  the  overestimation  introduced  in 
reaction  time  is  mitigated  by  the  square  root,  we  may  conclude  that  this 
overestimation  is  fairly  minor. 

A very  important  step  in  the  course  of  this  work  was  the 
identification  of  the  'final'  or  post-recombination  yields  in  terms  of 
the  'initial'  yields.  It  was  at  this  point  that  we  introduced  the 
approximation  of  globally  instantaneous  charge  recombination.  Although 
this  step  allowed  us  to  shrink  the  problem  down  to  more  manageable 
proportions,  greatly  economizing  the  computer  code,  as  was  noted  above. 

It  also  washes  out  a very  real  effect,  namely  the  effect  of  charge 
separation.  It  is  therefore  with  very  real  interest  that  we  look  forward 
to  the  upgrading  of  this  code  in  order  to  take  charge  separation  into 


account. 
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Later,  after  identifying  the  dependence  of  ion  reaction  time  scales 
on  medium  density  and  injection  rate— by  way  of  rather  simple 
considerations,  it  became  possible  to  quantitatively  identify  the 
dependence  on  medium  density  and  injection  rate  of  the  separation  of  ion 
and  neutral  time  scales  or  the  lack  of  it.  Specifically,  they  were 

found  to  merge  for  density  near  that  of  liquid  and  for  injection  rate  at 
about  the  early/late  chemistry  boundary. 

In  this  situation  charge  recombination  can  no  longer  be  viewed  as 
even  locally  instantaneous,  and  whether  there  is  charge  separation  or 
not,  it  becomes  necessary  to  accommodate  noninstantaneous  charge 
recombination  via  a split  time  step,  which  is  what  was  envisaged  in 
order  to  accommodate  charge  separation,  anyway.  The  fix-up  of  the  code 
is  the  same  in  either  case.  This  result  provides  all  the  more  impetus  to 
upgrade  DIFKIN  in  this  way. 

physical  relevance.  Without  physical  relevance,  the  diffusion 
kinetics  calculations  become  merely  a mathematical  exercise.  In  the 
interest  of  securing  physical  relevance,  the  heating  rate  by  the 
electron  beam  of  the  medium  has  been  estimated,  albeit  roughly,  without 
taking  into  account  energetics  of  the  reactions,  or  the  details  of  the 
energy  loss  process  in  the  initial  deposition. 

Experimental  verification.  This  is  an  appropriate  point  to  next 
turn  to  experimental  verification.  Certain  obstacles  seem  to  be  present 
here,  even  in  the  cases  in  which  overheating  does  not  appear  to  be  a 
problem.  These  are,  specif ically—limits  on  the  duration  of  a 
continuous  pulse,  limits  of  spatial  resolution  of  spectroscopic 
detectors  to,  perhaps,  the  0.1-1  mm  range,  possible  wind  effects  at  the 
tip  of  the  electron  gun  due  to  differential  pumping.  In  this  connection 
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the  technique  of  an  aluminum  window  may  be  of  use.  It  can  also  serve  as 
a downstepper  of  incident  energy  in  the  CSDA  limit  where  collimation  and 
monochromaticity  would  not  be  much  affected.  Boundary  effects  on 
spatial  yields  in  the  backplume  are,  of  course,  unavoidable  either  way. 

An  alternate  future  plan  of  work  has  been  proposed  which  may  prove 
useful  in  collaboration  with  Hanrahan  and  coworkers.  The  feasibility  of 
such  a plan  seems  within  the  realm  of  possibility. 

_^trapolation  to  liquid  density.  We  may  lastly  mention  one  aspect 
of  the  work,  the  extrapolation  to  liquid  density,  which  may  be  of  some 
significance.  In  the  first  place  these  cases  are  physically  realizable, 
since  heating  is  not  a problem  for  the  injection  rates  considered.  They 
have,  however,  two  deficiencies.  One  is  that  they  are  calculated  from 
spatial  yields  which  are  calculated  on  the  basis  of  gaseous  phase  energy 
loss  processes.  Referring  to  the  work  of  Kutcher  and  Green  (1976a),  we 
find  that  they  were  able  to  model  the  effect  of  the  plasmon  (the 
collective  excitation  mode  in  the  liquid)  on  energy  loss  in  liquid  water 
by  a broadening  of  the  general  oscillator  strength  peak  in  energy,  and 
found  a modest  effect  on  W value.  Referring  to  the  later  work  of  Turner 
_e^£l_.  (1982),  we  find  25-30%  differences  reported  in  calculations  of 
ratio  of  ionization  cross  section  to  total  inelastic  cross  section 
(ionization  plus  excitation).  A similar  difference  is  seen  in  the  peak 
of  mass  stopping  power;  however,  slowing  down  spectra  for  1 keV  and  10 
keV  electrons  are  not  seen  to  differ  greatly  over  several  orders  of 
magnitude.  One  is  led  to  the  conclusion  that  such  extrapolation  of 
spatial  yield  may  be  of  at  least  fair  accuracy.  This  first  deficiency  is 
perhaps  then  not  so  great. 
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Secondly,  the  aqueous  electron  (Kupperman,  1974;  Hanrahan,  1983) 
and  aqueous  0H“  and  their  reactions  (Kupperman,  1974)  have  not  been 
taken  into  account.  This  could  be  a more  serious  matter.  However,  it 
should  not  prove  difficult  to  encode  kinetics  for  a new  set  of  species 
and  reactions.  We  may  then  compare  with  results  obtained  using  only 
gaseous  phase  reactions.  Results  presented  here,  based  on  the  assumption 
of  local  instantaneity  in  the  case  of  the  kinetics,  and  global 
instantaneity  in  the  case  of  the  diffusion  kinetics,  have  introduced 
still  more  approximation.  As  was  pointed  out  above,  the  assumption  even 
of  local  instantaneity  in  the  case  of  liquid  water  is  on  the  very 
borderline  of  validity.  Given  these  differences,  it  would  seem  premature 
at  this  point  to  make  any  predictions  concerning  comparisons. 


Conclusion 

In  conclusion,  we  may  say  that  we  have  considered  in  rather 
thorough  detail  what  would  happen  if  a continuous,  steady  electron  beam, 
suddenly  turned  on,  were  injected  into  a target  of  water  vapor.  The 
time  scales,  principal  pathways,  the  effect  of  diffusion,  the  dependence 
on  medium  density.  Injection  rate,  temperature,  and  incident  energy  have 
all  been  considered  in  detail. 

One  thing  which  was  not  done,  but  which  would  have  provided  very 
valuable  information,  and  should  be  done  in  future  work,  would  have  been 
to  calculate  and  display  reaction  rates  and  reaction-specific  production 
rates  at  r=0  in  the  diffusion  kinetic  case,  for  comparison  with  the 
corresponding  pure  kinetic  rates.  This  information  can  still  be 
extracted  from  the  data  stored  on  tape  of  the  runs  which  were  made,  and 
should  prove  very  Interesting. 
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To  the  purpose  of  this  thesis  a computer  code  has  been  developed, 
which  is  capable  of  being  upgraded,  so  that  effects  so  far  not 
considered  may  also  be  investigated,  thereby  extending  the  applicability 
of  this  work.  The  fact  that  this  code  was  adapted  from  a vastly 
different  code  written  for  a vastly  different  purpose  (and  that,  along 
the  way,  this  original  code  had  been  adapted  in  a not  so  extreme  way  to 
a more  similar  problem)  attests  to  the  robustness  of  the  numerical 
machinery,  and  the  excellent  utilitarian  concept  in  the  design  of  this 
original  code. 

Because  the  code  is  time  dependent,  it  should  prove  equally  capable 
of  calculating  for  the  associated  problem  of  an  initial  pulse  beam. 

This  may  be  of  more  practical,  immediate  use  in  consideration  of  present 
experimental  capabilities.  Many  of  the  mathematical  tools  developed  and 
arrayed  here,  and  a great  deal  of  the  methodology  developed  here,  should 
prove  useful  in  such  a work. 


APPENDIX  A 

ADAPTATION  OF  THE  BOAT  CODE 

The  BOAT  code,  written  by  Sanford  Dash  and  Harold  Pergament  of 
Aeronautical  Research  Associates  of  Princeton  Associates,  Inc.  for  the 
U.  S.  Army  (Dash  and  Pergament,  1978)  to  calculate  turbulent  jet 
entrainment  with  chemistry  in  rocket  exhaust  plumes,  may  seem  like  an 
unlikely  place  to  start  in  order  to  base  a code  which  should  calculate 
diffusion  kinetics,  whether  time  dependent  or  steady  state,  for  electron 
beam  injection  in  a gas.  However,  this  code  does  many  things  and  can  be 
adapted  in  many  different  directions.  Finding  a meeting  place  half  way 
in  between  what  this  code  can  do  and  what  one  may  want  to  apply  it  to 
can  make  a difference.  In  order  to  better  appreciate  this,  the  place  to 
start  is  with  a brief  description  of  what  the  BOAT  code  is,  what  it 
does,  and  what  options  it  provides. 

In  its  most  glorious  form,  the  BOAT  code  consists  of  three  job 
steps:  PROCES,  SCIPPY  and  BOAT.  The  first  step  reads  in  data  and  puts 

together  a chemical  system  of  species  and  reactions  as  directed, 
referencing  its  list  of  species  and  their  thermodynamic  data  (for  the 
purpose  of  calculating  enthalpies,  specific  heats  and  Gibbs  free 
energies  contained  in  the  form  of  polynomial  fits  based  on  JANNAF  Table 
data)  and  its  list  of  reactions  and  reaction  rate  parameters. 

The  second  step  calculates  supersonic  aerodynamic  inputs  for  the 
third  step,  if  applicable,  via  what  is  known  as  an  "overlaid" 
procedure,  and  can  be  bypassed.  For  our  purposes  then,  it  can  be 
bypassed.  The  third  step,  which  we  are  most  interested  in,  calculates 
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turbulent  jet  entrainment  (modelling  by  an  effective  diffusion)  with 
chemistry  in  the  flowfield  of  the  exhaust  plume.  It  can  therefore  be 
said  to  calculate  steady  state  diffusion  kinetics  in  an  axial  flow 
field. 

The  first  point  of  arrival  to  which  it  was  desired  to  take  this 
code  was  to  apply  it  to  the  problem  of  axisymmetric  laminar  flames,  both 
diffusion  and  pre-mixed  (see  Kanury,  1975). 

The  next  requirement  was  to  debug  the  code  of  a few  lingering 
errors  and  typos,  of  which  two  had  to  do  with  the  encoding  of  two  of  the 
five  reaction  type  pairs  which  are  encoded  in  BOAT: 

(1)  A + B -»■  C + D 

(2)  A+B  + M -»■  C + M 

(3) A  + B -»■  C + D + E 

(4)  A + B C 

(5) A  + M C + D + M 

where  (6)  - (10)  are  the  corresponding  two-way  reactions.  These  four 
reaction  types  (3,8;  5,10)  are  not  represented  in  the  six  standard 
chemical  systems  built  into  the  code,  which  leads  me  to  believe  that 
they  had  never  been  checked. 

The  next  and  final  step  was  to  adapt  the  code  to  this  particular 
problem.  The  numerical  machinery  of  the  BOAT  code  (job  step  three) 
accomplishes  the  following.  The  hydrodynamic  equations  of  axial  flow 
combined  with  radial  diffusion  in  streamline  coordinates  govern  the 
conservation  of  angular  momentum,  energy  or  enthalpy,  and  species.  If 
nonlaminar  viscosity  is  opted  for,  then  equations  governing  conservation 
of  turbulent  kinetic  energy  and  turbulent  dissipation  are  also 
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included.  However,  in  the  laminar  mode  these  equations  are  not 
included.  The  three  equations  of  interest,  as  they  appear  in  BOAT,  are 

— = - — + — — Ta^'i  . ^ 
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where  u is  axial  velocity,  v is  radial  velocity,  x is  axial  coordinate, 
is  streamline  coordinate  or  streamline  function,  T is  temperature, 
y is  coefficient  of  viscosity,  p is  density,  p^  is  ambient  edge 
density,  and  g is  acceleration  due  to  gravity,  Cp  is  specific  heat,  F.  = 
X^/W  where  is  mole  fraction  of  i-th  species  and  W is  mixture 
molecular  weight,  p is  pressure,  h^  is  enthalpy  of  i-th  species,  r is 
cylindrical  radial  coordinate,  w^  is  net  rate  of  production  of  the  i-th 
species,  and  Le  is  Lewis  number  and  Pr  is  Prandtl  number  (see  Kanury, 
1975),  and  the  quantity  A is  given  by 


A = 


y 


2 

pur 

’1' 
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The  tranf ormation  from  (x,r)  to  (x,i|;)  is  given  by 


dip 

<1^  -^  = pur 


(A. 5) 
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I ^ 

(A.6) 

In  the  case  of  constant  pressure  mixing  the  pressure  gradient  term  in 
the  axial  momentum  and  energy  equations  drops  out. 

The  mystery  term  Z is  stated  by  the  authors  to  be  (9u/9tj))^,  but 

its  stated  discretization  corresponds  to  9^u/9i|;^.  It  is  encoded  as 
2 

u . By  comparison  with  other  standard  treatments  (for  example,  see 
Miller  and  Kee,  1977)  it  is  apparent  that  this  term  represents  a small 
correction  to  stagnation"  enthalpy  and  can  be  neglected  at  lower 
velocities. 

This  set  of  equations  is  parabolic  in  the  axial  coordinate,  and 
therefore  a parabolic  marching  technique  is  implemented  to  solve  it  for 
given  nozzle  conditions  ('initial'  conditions)  and  given  mixing  layer 
boundary  conditions,  i.e.  at  the  edge  which  moves  outward  radially  with 
increasing  x.  The  radial  coordinate  is  transformed  via  the  von  Mises 
transformation  to  the  streamline  coordinate,  ip,  and  the  streamlines  are 
evolved,  via  enthalpy  conservation,  so  that  the  radial  grid  is  made  to 
follow  along  constant  streamlines,  equally  spaced  in  streamline 
function,  \p.  The  streamline  evolution  is  accomplished  via  grid-point 
redistribution  technique  which  incorporates  a cubic  interpolation 
rule.  This  grid-point  redistribution  technique  is  one  of  the  state  of 
the  art  aspects  of  the  code.  However,  its  usefulness  is  limited  to  the 
flowfield  expansion  context. 

We  may  note  in  passing  that  as  a result  of  transformation  to 
streamline  coordinates,  the  equation  of  continuity. 
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9pu  1 3 . . 

-^  + 7^(prv)  = 0 (a.7) 

is  solved  for  automatically. 

The  contribution  of  the  mass  flux  term  in  the  energy  equation  is 
taken  into  account  via  a predictor-corrector  method  in  enthalpy, 

inverting  temperature  from  enthalpy  at  each  step,  based  on  the  previous 
value. 

The  kinetic  production-destruction  terms  are  calculated  via  an 
implicit  linearized  approach.  This  technique  is  known  to  guarantee 
numerical  stability.  The  implicit  linearization  is  given  by 


- (rpj  + - (F.)J 
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Here  the  subscripts  indicate  axial  station.  Since  the  values  at  station 
n are  known,  the  values  at  station  n+1  are  to  be  solved  for.  They  are 
therefore  said  to  be  defined  implicitly.  They  are  obtained  by  solution 
of  the  simultaneous  system  of  linear  equations  which  the  linearization 
technique  implies. 
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Application  to  Laminar  Flames 

The  adaptation  of  the  BOAT  code  to  the  laminar  flame  problem 
required  the  following  modification  of  the  BOAT  code  numerical 
machinery.  In  place  of  constant  values  of  Lewis  and  Prandtl  numbers, 
temperature  dependent  values  were  required.  The  values  implemented  in 
the  BOAT  code  for  these  dimensionless  numbers  were  constants  of  order 
unity.  Since  they  are  temperature  dependent  and  inside  the  streamline 
derivative,  the  finite  spatial  differences  had  to  be  rewritten  to 
include  their  spatial  variation  (as  a result  of  nonzero  temperature 
gradient) . 

Noting  the  standard  definitions  of  these  dimensionless  numbers 
(Kanury,  1975),  it  becomes  clear  that  the  ratio,  Le/Pr,  appearing  in  the 
mass  flux  term  in  equations  (A. 2)  and  (A. 3)  is  incorrect,  and  should  be 
replaced  by  1/PrLe.  It  is  worth  going  over  the  details.  The  standard 
definitions,  taken  from  Kanury  (1975),  are 

Pr  5 v/a  = (u/p)/U/Cp)  (A.io) 

Le  = a/D  = (X/Cp)/D  (A.n) 

The  viscosity  factor,  A = ppur^/i()  = vp^ur^/i|),  and  the  thermal  factor, 

2 

CpA/Pr  = Xpur  /ijj,  are  correct  as  they  stand.  The  correct  mass  flux  or 
diffusion  factor  is  A/PrLe  = Dp^ur^/iJ;,  not  LeA/Pr. 

Finally,  it  should  be  noted  that  in  the  laminar  flame  application 
constant  pressure  mixing  is  assumed,  as  a result  of  which  the  pressure 
gradient  terms  in  Eqs.  (A.l)  and  (A. 2)  drop  out. 
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This  modified  BOAT  code  was  implemented  to  solve  for  various  flame 
conditions:  both  H2~air  and  CH^-air,  both  diffusion  and  pre-mixed. 
Interesting  results  were  obtained,  especially  in  comparison  with 
published  results  of  Miller  and  Kee  (1977),  Kee  and  Miller  (1978),  and 
with  our  own  experimental  data  (Pamldimukkala  et  al.,  1983).  However, 
the  author's  involvement  on  this  project  was  terminated  for  lack  of 
funds.  Fortunately  he  could  be  supported  by  the  radiological  grant  by 
applying  this  code  to  the  radiological  problem. 


Application  to  the  Radiological  Problem 

The  following  problem  is  now  posed.  Adapt  the  BOAT  code  to  the 
problem  of  diffusion  kinetics  for  an  electron  beam  incident  in  water 
vapor.  This  problem  is  parabolic  in  time  coordinate,  and  elliptical  in 
spatial  coordinate.  The  BOAT  code,  which  in  effect,  calculates  steady 
state  diffusion  kinetics  in  an  axial  flowfield,  solves  a parabolic 
problem  in  two  dimensions. 

The  general  radiological  problem,  being  axlsymmetric,  involves  one 
parabolic  dimension,  time,  and  two  elliptical  dimensions,  namely,  the 
two  spatial  coordinates,  say,  cylindrical  coordinates. 

One  logical  choice  therefore  greets  us:  either  pose  the 

radiological  problem  in  one  (radial)  spatial  coordinate,  or  expand  the 
BOAT  code  to  cover  a two  dimensional  half-azimuthal  grid.  The  decision 
was  made  to  go  the  first  route. 

The  idea  occurred  of  constructing  a generalized  orthogonal 
coordinate  system  of  radial  type,  ellipsoidal  in  shape,  in  order  to 
allow  calculation  in  a somewhat  more  faithful  spatial  form.  The 
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author's  idea  was  that  if  the  contours  were  made  to  be  spherically 
asymptotic,  therefore  realistically  reflecting  the  trend  of  diffusion, 
in  fact,  the  trend  of  the  initial  radiological  deposition  process  as 
well,  that  by  neglecting  tangential  diffusion  a one  dimensional 
approximation  to  two-dimensional  diffusion  would  be  possible. 

Such  contours  were  actually  encoded,  and  this  work  is  described  in 
Appendix  B.  However,  two  things  happened:  (1)  an  esoteric  impasse 
developed  in  the  checking  out  of  the  coordinate  system.  Calculated  pure 
diffusion  approach  to  steady  state  showed  tolerable  drift  at  the  origin 
for  the  polar  coordinate  ranging  from  the  radial  direction  to 
infinitesimally  close  to  the  axial  direction  (n  = 0.999999).  However, 
in  going  over  to  polar  contour  along  the  axis  (n  = 1),  a jump  occurred 
which  the  author  could  not  explain  (see  Fig.  A.l).  (2)  Typical 

deposition  contours  at  incident  energies  up  to  and  around  1 keV  proved 

to  be  of  sufficiently  low  eccentricity  so  that  a spherical  approximation 
would  suffice. 

It  was  then  decided  to  go  ahead  with  this  work  in  the  spherical 
approximation.  It  should  be  noted  in  passing  that  it  should  also  be 
easy  to  adapt  the  code  to  the  case  of  cylindrical  geometry.  As  to  the 
general  case,  it  would  seem  that  after  the  esoteric  impasse  is  overcome, 
the  generalized  radial  coordinate  system  could  prove  useful  in  the 
implementation  of  a two  dimensional  grid.  The  advantage  that  it  could 
provide  would  be  the  following.  With  contours  aligned  with  the  source 
in  the  core  region  and  asymptotically  spherical  in  the  limb,  an  optimal 
representation  would  be  provided  at  minimal  expense,  i.  e.  with  a 
minimal  polar  grid. 
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Figure  A. 1 . Pure  diffusion  time  course  at  ^ = 0 for  n = 0,  .999999,  and  1. 

Results  of  numerical  code  Implementing  one  dimensional  approximation. 


283 


Subroutine  Organization  of  DIFKIN 

In  this  section  the  subroutine  structure  and  function  of  DIFKIN  are 
summarized  in  order  to  provide  a user's  guide  to  the  code  and  to 
document  the  adaptation  from  BOAT.  With  the  first  aim  uppermost  in  the 
author's  mind,  the  summaries  here  have  been  made  as  complete  and  as 
clear  as  possible. 


Job  Step  1: 


MAIN  reads  in  input  data:  chemical  system  index  numbers, 
physical  input  data,  control  parameters  (pointers  and 
quantifiers),  then  writes  the  chemical  system  along  with 
regurgitation  of  other  inputs  to  temporary  disk  file  which  is 
read  by  subroutine  SA  in  job  step  2.  Calls  SYSMOD  (optionally 
including  entry  point  RRMOD)  and  DABANK  (including  entry  points 


LOCOB  and  REACT). 

2.  SYSMOD  constructs  chemical  system  (references  BLOCK  DATA, 
entry  point  RRMOD  rescales  chosen  reaction  rates. 

3.  DABANK  — reads  in  species  ID's,  mol.  wts.,  enthalpies  of 


formation,  and  fitting  parameters  for  thermodynamic  functions 
from  card  image  file  "CHEMRAD"  (named  by  author).  It  then  reads 
in  reactions,  i.  e.  species  entries,  reaction  type,  and  reaction 


rate  coefficient  parameters. 


Job  Step  2: 

1.  MAIN  — "home":  1.  calls  SA  which  reads  in  data  from  temporary 

disk  file  written  to  by  MAIN  (job  step  1),  2.  calls  PREP  (taken 
bodily  from  Rio's  code  DIFKT  (Rio,  Ph.  D.  thesis,  1983))  which 
reads  in  input  data  for  calculation  of  spatial  yields,  namely 
plume  integrals  and  spatial  yield  spectra  parameters,  input  3. 


loops  over.  4.  calls  SB,  the  Initialization  subroutine.  (SB 
initializes  the  grid,  the  concentrations,  and  the  yields),  5. 
loops  on  Ml,  the  marching  step  driver,  until  it  is  flagged.  It 
then  calls  01  (entry  point  in  OT),  which  writes  permanent 
restart  disk  file,  and  then  goes  on  to  next  point,  if  the 
last,  then  stops. 

system  of  simultaneous  linear  equations  solver.  Gauss- 
Jordan  pivot. 

IC  — assigns  initial  conditions:  value  of  1.  for  H2O,  value  of 
10  for  all  others  (concentrations  in  units  of  the  H2O 
concentration,  mole  fractions  if  dilute  regime). 

PREP  — described  in  MAIN. 

YI  — extracted  from  Rio’s  code,  DIFKT.  Calculates  spatial 
species  yields  via  branching  ratios  over  state  yields,  folded 
into  spatial  plumes  of  spatial  yield  spectrum  (see  Eq.  2.25). 
Post-recombination  reassignment  of  spatial  species  yields  is 
then  made  (this  thesis,  see  Eqs.  4.3a-g,  4.5  and  4.6). 

NR  — "Newton-Raphson,  Euler".  Calculates  5"grid.  Depends  on 
value  of  q. 

SIMSON  ~ Simpson  rule  interpolation.  Two  versions:  previous 

dangerous  version  with  convoluted  search  logic,  later  safer 
version  with  bubble  search. 

CC  calculates  entries  in  column  vector  and  matrix  implied  by 
linearization  of  kinetic  production-destruction  terms.  The 
encoding  of  this  subroutine  is  extremely  compact. 


This  is 
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advantageous  to  computer  efficiency  but  renders  its  content 
totally  opaque  to  the  user  until  after  several  hours  of 
combinatorial  and  algebraic  deciphering. 

9.  GX  calculates  Gibbs  free  energy  for  each  species. 

10.  EF  — (function  subprogram).  Calculates  Arrhenius  factor  for 

backward  reaction  rate  if  two  way  reaction  type  is  opted,  calls 
GX. 

a.  calculates  forward  (if  one  way)  or  forward  and 
backward  (if  two  way)  kinetic  production  and  destruction  terms, 
then  calls  EF  if  two  way,  then  calls  CC  which  uses  the  above 
terms  as  inputs  (see  above);  b.  before  returning,  calls. 

12.  S3  Integrator  engine:  calculates  finite  difference  diffusion 

Laplacian  and  implicit  chemistry  via  call  to  SL  using  Inputs 
from  CC.  Evolves  all  points  except  edge  point.  Center  point 
(i.e.,  origin)  is  evoled  separately  via  third  difference 
extrapolation. 

13.  Ml  marching  step  driver:  1.  calculates  time  step:  only 

diffusion  unless  chemistry  override  is  specified,  if  that  is, 
then  takes  minimum  of  the  two  until  diffusion  time  steps  crosses 
plateau  value  of  chemistry  time  step,  after  which  for  the  sake 
of  economy  it  takes  the  geometric  average  (calls  function 
subprogram  CHEM  to  evaluate  chemistry  time  step);  if  pure 
kinetic  run  time  step  equals  chemistry  time  step,  2.  decides 
whether  to  print  out  results  at  this  time  station  or  not  on 
basis  of  either  (a)  discrete  blocks,  'unlimited'  number  (choice 
in  original  BOAT  code  was  binary),  or  (b)  newly  coded 

pj-  “ ^pr^’  IN  to  output  updated 


exponential  rate:  At 
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geometry  of  the  grid,  regurgitate  miscellaneous  physical  input 
data,  and  output  the  spatial  yields  at  the  present  grid  points 
(if  pure  kinetic  run,  skip  all  this  after  first  time  station), 

4.  edge  point  extrapolation  (second  difference  in  logarithm), 

5.  calis  M2,  and  6.  checks  for  final  time  exceeded,  if  so 
returns  to  loop  in  MAIN. 

CHEM  (function  subprogram):  calculates  "chemistry  friendly" 

times  step.  See  (Eq.  5.1). 

M2  two  functions:  1.  check  for  negative  "mole  fractions" 

(precisely,  below  -10  ^*^).  if  so,  halve  the  time  step  and  try 
again,  if  time  step  falls  below  minimum,  print  warning  and 
stop,  if  "mole  fraction"  exceeds  -lO'^^  but  is  exceeded  by 
10  truncate  to  10  and  carry  on,  2.  grid  point 
redistribution:  evolve  5^,  redistribute  grid  in  discrete  steps 

depending  on  whether  or  not 


new  5 (t) 



value  at  last  print  out 


> c 

int 


then  call  SIMSON  to  Interpolate  in  the  logarithm  of  the 

concentration  over  the  transformed  coordinate, 

C = ln(C/C  + 1). 

o 

IN  ~ writes  out  header  and  inputs.  Adapted  to  write  out 
updated  geometric  information,  misc.  physical  inputs,  and  yields 
at  updated  grid  points  when  called  from  Ml  (see  Ml). 

reads  inputs  from  temporary  disk  file  written  by  MAIN 
in  job  step  1,  2.  reads  in  diffusion  coefficients  at  300  K and  1 
atm  from  unit  5,  calculates  diffusion  coefficients 


at  specified 
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temperature  and  density,  3.  assigns  numerical  indices  to  each 
reaction  to  indicate  species  and  positions  (this  aspect  of 
original  BOAT  code  required  some  debugging). 

18.  SB  ~ initializes  (1)  (2)  Gibbs  free  energies,  (3)  certain 

control  parameters  which  change  value  through  cycles,  (4) 
initializes  grid  (calls  NR),  concentrations  (calls  IC)  and 
yields  (calls  YI),  and  (5)  converts  (overwrites)  initial  "mole 
fractions  (molar  concentration/molar  concentration  of  H2O)  to 
working  units  (adapted  term  of  Dash  and  Pergament),  namely, 
mole  fraction  /mol.wt.  of  H2O.  "Working  units"  in  original 
code  were  true  mole  fractions/mixing  mol.  wt., 

L9.  OT  'output'  writes  out  concentrations  as  such  or  scaled 

by  H2O  concentration  [and  other  Information  to  printer  and  to 
disk  file],  entry  point  01:  writes  out  RESTART  file  (disk), 

entry  point  02  (called  from  S2  if  opted  for):  write  out 

reaction  rates. 


Grid  Point  Redistribution  and  Time  Step 
The  grid  point  redistribution  scheme  existing  in  the  BOAT  code  was 
found  to  be  inappropriate  to  our  purposes.  A new  one  was  therefore 
written.  Although  a cubic  scheme  is  preferable  to  handle  inflections, 
if  a transformed  variable  can  be  found  in  which  the  concentrations  are 
inflection  free  or  very  close  to  it,  then  a Simpson  rule  is 
sufficient.  A Simpson  rule  Interpolating  subroutine  was  written. 

The  distribution  of  the  radial  coordinate  (?)  grid  is  our  next 
consideration.  In  the  interest  of  optimal  bin 


averaging  we  require  that 
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the  packing  of  the  grid  reflect  the  evolving  concentration  profiles.  In 

order  to  accomplish  this  we  set  up  an  equally  spaced  grid  In  a 

transformed  coordinate,  ? = ln(?/£  (t)  + 1)  where  5 (t)  Is  a monotonlc 

o o 

Increasing  function  of  t through  In  t,  chosen  empirically,  the  purpose 
of  which  Is  to  control  the  evolution  In  time  of  the  grid  redistribution, 
and  through  It,  the  time  dependence  of  the  time  step  (see  Eq.  A.l),  As 
far  as  only  the  diffusion  Is  concerned,  the  trade~off  Involved  here  Is 
between  cost  and  the  accuracy  of  discrete  covering  of  the  central 
core.  Accuracy  In  this  Is  not  critical  to  the  overall  accuracy  of  the 
calculation,  and  may  be  compromised  to  some  extent,  with  negligible 
Impact  upon  overall  accuracy. 

The  grid  endpoints  In  ^ are  fixed  at  0 and  an  edge  point, 

= 36/8d  (where  d Is  medium  density).  The  evolution  at  the  edge 
point  Is  accomplished  through  second  difference  extrapolation  In  the 
logarithm.  The  value  = 36/ 8d  Is  chosen  since  at  this  point  the 
steady  state  diffusion  profile  has  fallen  off  to  5%  of  Its  peak  value  at 
the  center,  for  y(r)  = e 

The  discretization  of  the  Laplaclan  at  the  origin  Is  accomplished 
through  a third  difference  extrapolation  and  half-bin  assignment  (In 
order  not  to  divide  by  zero). 

Discretization  of  the  derivatives  with  respect  to  ^ Is 
accomplished  through  transformation  of  the  derivatives  with  respect  to 
the  equally  spaced  coordinate, 

9u  _ 

35  ~ d5  3?  (A. 10) 

2 2 

y.  _ d^  rd /'d5'\3u  ^ 

3^2  d?  ''d?^  ^ >^d5J 


(A. 11) 
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Since 


rii'i  _ dc  _ 1 

dg  ■ dC  ■ " 5 + 5 (A. 12) 


for  the  transformation,  ^ - ln(g/g^  + 1),  the  transformation  formula 
for  the  second  derivative  is  given  by 


(I  + 


r 8 U 8u-i 

^,2  ■ acJ 


(A.13) 


wherefore  the  interpolated  second  difference,  is  given  by 


A^u.  = 

X 


- 2 +y^)u.  ,1. 

(g^  + 5^)^  ^ 1 2 1-1 J 


(A.14) 


The  step  size  determination  already  existing  in  the  BOAT  code 
embodies  the  standard  stability  criterion  (see  Kupperman  and  Belford, 
1961a), 


At  = 4 (h, A5)^,  /D 

J i min  max 


(A.15) 


expressing  diffusion  control,  with  repeated  halving  of  this  value  if 
negative  concentrations  result.  The  generalized  difference  h^^A^ 
(explained  in  Appendix  B)  reduces  to  Ar  in  the  spherical  case. 

A chemistry  friendly  override  time  step  has  been  added  to 
accommodate  the  chemistry  for  pure  kinetic  runs  and  diffusion  kinetic 

runs  in  the  early  or  borderline  chemistry  regime. 

« 

In  kinetic  codes  based  on  the  Gear  method  (Hindmarsh,  1974),  the 
time  step  is  determined  on  the  basis  of  the  derivatives  of  the 
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solution.  We  find  that  a simple,  empirical  approach,  proves  adequate  in 
the  context  of  the  calculation  of  kinetic  production  - destruction  terms 
via  implicit  linearization  (Eqns.  (A. 8)  and  (A.9))  employed  in  the  BOAT 
code.  Specifically  we  find  that  the  concentration  time  courses  became 
stabilized  (i.e.  converged  in  function  space)  for  At(t)  everywhere 
less  than  or  equal  to 


At 


stbl 


0. It,  t < T 

C 

0.  1 T , t > T 

C C 


(A. 16) 


where  is  reaction  time  (see  Chapter  III,  section  2).  It  was  also 

found  that  this  requirement  could  be  relaxed  somewhat  for  t >>  r with 

c 

slow  logarithmic  growth.  This  was  very  advantageous  for  economy  in  the 
context  of  the  full  diffusion  kinetic  calculations.  These  results  are 
illustrated  in  the  selection  of  kinetic  time  plots  presented  in  Fig. 


A. 2. 
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Figure  A.  2.  Kinetic  time  courses  for  various  choices  of  At(t). 
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Figure  A. 2 


continued 


APPENDIX  B 


GENERALIZED  ORTHOGONAL  COORDINATE  SYSTEM  OF  RADIAL  TYPE 
An  attempt  was  made  to  encode  a generalized  radial  grid,  and  in 
fact  came  within  a very  small  fraction  of  success.  However,  it  became 
apparent  that,  aside  from  an  esoteric  problem  that  arose  (see  discussion 
in  Appendix  A),  such  a way  to  go  would  only  really  pay  off  if  it  were 
carried  all  the  way  to  a two  dimensional  grid.  Where  it  was  left  at  was 
not  very  far  from  there.  For  the  sake  of  documentation  and  reference, 
this  work  is  presented  here. 

Let  the  orthogonal  coordinates  be  5 and  p,  where  ? is  the 
generalized  radial  coordinate  and  p is  the  generalized  polar 
coordinate. 

first  choose  to  consider  the  case  of  the  isocontours  being 
concentric  prolate  ellipses: 

5 = / + (z  - z^)^/b^(C)  . (B.l) 

The  axis  ratio  parameter  b(5)  is  chosen  to  be  a smooth,  slowly  varying 
function,  monotonically  decreasing  from  a maximum,  not  necessarily  a 
stationary  maximum,  at  5 = 0 as  ^ Referring  to  Appendix  C where 

steady  state  axlsymmetric  diffusion  is  calculated,  a lesson  is  gained 
concerning  the  limits  of  ‘intuition":  the  maximum  of 

b(C)  at  § = 0 for  the  steady  state  solution  for  source  term  of  the  form 


9QT 
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y(c.z)  . Ae-’' 

(where  is  a constant)  is  not  stationary. 

It  may  be  noted  that  if  5 is  thought  of  as  a function  of  p and 
z,  it  may  also  be  thought  of  as  a functional  of  b.  We  may  therefore 
write  5 = ?[b](p,z). 

The  perpendicular  slope  follows  at  once  from  Eq.  (B.l): 

rlP)  = b^(g)P 

^9z'^5  I z— z r * (B.2) 

We  now  introduce  the  adiabatic  approximation,  where  we  consider  b(^)  to 
be  locally  constant.  Equation  (B.2)  then  locally  integrates  to 

k(5,n)  = |z  - z^l*’  ^^Vp  . (B.3) 

We  now  introduce  the  local  coordinate,  p(5,q)  which  may  be  taken  as  any 
suitable  monotonic  extended  analytic  function  of  k.  In  order  to  regain 
the  conventional  definition  in  the  spherical  limit,  n = cos 9,  we  set 
q = k//l  + k for  p > 0,  and  q = ±1  for  p = 0,  z 0.  We  now  define 
the  global  or  governing  orthogonal  coordinate,  q,  for  the  sake  of 
definiteness,  by  the  boundary  condition,  q(0,q)  = q.  (The  other 
suitable  definition  would  be  q(«.,q)  = q.)  it  can  be  shown  that 
|q|  is  monotonlcally  decreasing  in  with  maximum  |q|  = |q(0.q)| 

and  minimum  |q|^.^  = |q(»,q)|  for  0 < |q|  < 1,  and  stationary  at  the 
turning  points,  |q|  = 1,  0.  We  may  now  construct  the  q-contours. 

This  is  done  by  Integrating  the  first  order  normal  differential  equation 
given  by  Eq.  (B.2)  where  the  dependence  of  the  right  hand  side  on  p and 
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z Is  implicitly  defined  with  the  help  of  Eq.  (B.l).  This  is  easily  done 
numerically  by  the  Euler  method,  by  which  we  arrive  at  a span  for  p and 
5 over  an  interpolatory  grid  in  z for  a given  value  of  p.  Contours 
of  C and  p are  shown  for  b(0  given  by 

b(?)  = 1 + 6 e (B.4) 

where  6 = 0.6  and  = 7.4  in  Fig.  B.l  where  it  is  also  used  to  display 
approximate  characterization  of  the  steady  state  solution  for  pure 
diffusion  for  source  term  with  b = 2. 

The  Laplacian  operator  in  this  generalized  coordinate  system, 
following  Morse  and  Feshbach  (1953),  is  given  by 


V^u 


1 3 fV  Ju>, 

h^h^p  85  ^h^  85^ 


1 J_ 


3?  '■■'2 


8p 


(B.5) 


where  h^  and  h2  are  the  scale  factors  of  5 and  p.  Following  Morse 
and  Feshbach,  we  note  that  they  are  given  by 


■i  ■ ^ - (H) 


95' 


>2  - ■'  (■^)"  + (^) 
8p  8p 


8z^2 


(B.6) 


The  adiabatic  approximation  is  invoked  in  the  evaluation  of  the  partial 
derivatives  by  implicit  differentiation  of  Eqs.  (B.l)  and  (B.3): 


9p/95  = 5p/5^[b^]  8z/85  = 5(z-z^)/5^ [b^  ] 


(B.7-8) 


z-z 


2-b 


b^5^[b^] 


3 I |l-b^  3 
..  P I z-z  I s 

9z/8p  = — — — 

b^5  [b^] 


8p/8p  = - 


(B.9-10) 
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Figure  B.l.  5 and  n contours  for  b(5)  given  in  Eq.  B.4  with 
<S  = 0.6  and  = 7.4,  ^ contours  are  isocontours 
of  approximate  solution  to  steady  state  diffusion 
equation  for  = 2 in  source  term  given  by  Eq.  C.9. 
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where  b - b(0,  and,  observing  that  n//l  - = k(?,n),  it  foiiows  that 

2 

S = dk/dn  = (k/n)^  = (|z  - /pn)^.  it  follows  that  h^  = c/?[b^] 

and  h^  = |z  - p^s^/b^?[b^]. 

In  order  to  encode  Eq.  (B.4)  the  special  cases  5=0  and  |n|  =1 
or  0 for  5 > 0 must  be  considered  separately  because  of  L'Hopital 
limits.  The  limiting  form  of  the  Laplacian  is 


,2^  . a,  ^ + ^2  |u 
1 3j2  5 35 


(B.ll) 


(in  the  5 term)  where  a^  = 1 and  a^  = 1 + 1/b^  for  p = 0,  and 
2 

= 1/b  and  = 2 for  |n|  = 1,  5 > 0 and  for  p ^ 0,  5 = 0. 

The  sense  of  the  values  of  a^  and  a2  is  made  clear  by  considering 
the  extreme  possible  limiting  cases  of  the  geometry  of  the  problem: 
b -»■  0 and  b -»■  «.  In  the  limit  b ->■  0 the  metric  becomes  locally  flat 
in  the  axial  direction  in  the  neighborhood  of  the  axis,  and  therefore, 
must  equal  1/b  for  scaling  and  dominate  a2  for 
limit  b =0  the  metric  becomes  locally  cylindrical 
of  the  z=0  plane,  and  therefore  in  this  limit  a^  = 
n = 0.  In  the  case  b = 1 the  spherical  metric  is 


I h I = 1 . In  the 
in  the  neighborhood 
a.2  = 1 for 
regained. 


APPENDIX  C 

AXISYMMETRIC  DIFFUSION— STEADY  STATE 
In  order  to  justify  spherical  approximation  it  is  useful  to 
calculate  pure  diffusion  in  the  general  axlsymmetrlc  case.  In  the 
general  axisymmetric  case 


R - / r - 2rr^[pp^  + / (l-y^)(l-p^^)cos(l)^  ] + r^^  (C.l) 

where  p = cos0  and  = cos6^,  the  Green's  function  integral  (Eq.  3.24) 
is  analytic  in  <[i^,  and  reduces  to 


c(r, 9,t) 


t dt 
' o 

0 (2/ttDt)^ 


X 


-2rr  cos(0-0 
o o 


) 


+ r ^)/4Dt 

° (1  + F(p)) 


(C.2) 


where  F(p)  = e ^I^(p)  - 1 and 


P = (2rr  /4Dt)/  (l-p2)(i-p 
^ o 


Iq(p)  Is  the  spherical  Bessel  function  of  order  0.  The  time  dependent 
problem  reduces  to  a triple  Integral  in  r^,  p^  and  t^.  Although 
numerical  calculation  of  the  contribution  to  this  integral  in  the  region 
t^  not  too  close  to  t is  straightforward,  it  becomes  tricky  in  the 
region  too  close  to  t because  the  singularity  of  the  Green's  function 


in  T. 
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Rsther  than  follow  this  out,  we  turned  to  the  steady  state 
solution,  which  is  more  easily  evaluated,  and  provides  suitable 
information  about  the  broadening  effect  of  diffusion  on  the  deposition 
contours  laid  down  by  the  source. 

Expressing  the  diffusion  Laplacian  in  spherical  coordinates,  the 
diffusion  equation  becomes 


i_  (',.2  ic-j  j.  i_  ^ r.,  2.  8c>|  y(r,u) 

2 8r  9r''  ^ ^2  3y  ^ ’ aJ  + n " 


d]i'  ' D 


(C.3) 


Expanding  in  Legendre  polynomials  in  p (Jackson,  1967),  we  set 
“ 00 
c = c(r,p)  = E c^(r)P^(p)  and  y = y(r,p)  = E y (r)P  (y)  where 


1=0 

y.(r)  is  given  by 


£=0 


Yg^(r)  = / ^y(r,u')P^(y’)dp'. 


(C.4) 


Noting  the  identity 


d 9 

^[(1-y  ) = 0, 


(C.5) 


it  follows  from  completeness  that 


9c , 


i--L_rr-2_A^  Uj^+1)  „ 

2 9r^  9r  2 *"Jl  D 0,  £ - 0,1,2, 


(C.6) 


Expanding  the  derivatives.  Equation  (C.6)  simplifies  to 


c"  + — c' ^.(  ^+1 ) Z 

Z r ^Z  2 ^Z  D 

r 


(C.7) 
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The  general  solutions  of  the  associated  homogeneous  equation  are 

^ _ j^_l 

= r and  = r . Noting  the  boundary  condition  at  infinity, 
using  the  method  of  variation  of  parameters  (Birkhoff  and  Rota,  1961), 
we  obtain 


c/r)  - 


Jl+1 


1+2 

r P y.(p)dp  oo  y (p)dp 

/ -rd ^ ^ 


0 


r (2£+l)p 


H-1 


(C.8) 


The  problem  thus  reduces  to  a single  integral  per  r value  plus  a basis 
set  of  single  integrals  over  p per  r value.  Values  of  c(r,y)  at  any 
y then  follow. 

We  have  calculated  c(r,y)  for  the  case  of  a concentric  elliptical 
source. 


o / ^ 2 2 

y(r.„)  . •* 

(C.9) 

(e  = / 1 - 1/b 

o 


The  contours  of  the  solution  maintain  the  shape  of  the  source 
(essentially)  near  the  center,  and  gradually  asymptote  towards  spherical 
at  large  distance. 

This  behavior  may  be  conveniently  represented  to  a high  degree  of 
accuracy  by  representing  c(r,y)  by  the  form 


c(r,y)  = c(5) 


(C.IO) 


where  the  generalized  radial  coordinate,  is  defined  implicitly  by 


the  relation 
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C = / + (z  - z^)^/b^(5)  (C.ll) 

2 

where  p = r/  1 - y and  z = + ry.  The  function,  b(C),  Is  the 

axis  ratio  parameter  of  the  half-elliptical  contour  in  the  p-z 
aximuthal  half-plane  defined  by  5 “ constant. 

The  5 dependence  and  the  trend  in  of  b(^)  is  interesting*  The 
it  approacheS“and  apparently  reaches"at  the  center  is  less 
than  b^  and  grows  sublinearly  in  b^-1.  If  this  maximum  is  stationary, 
it  is  so  only  within  a very  small  ball.  The  5 dependence  of  b(5)  is 
very  well  represented  by  the  exponential  form 

b(C)  = 1 + 6e  . (C.12) 

The  dependences  of  6 and  on  b^  are  adequately  represented  by  the 
forms  6 = 6^(b^  - 1)P  and  + 5^^(b^  - 1)  where  6^  = .59,  p = 

^bo  ^ ^bl  1*12.  Contour  graphs  of  the  approximate 

solution  to  the  steady  state  diffusion  equation  for  b = 2 are  plotted  in 
Fig.  (B.l)  along  with  calculated  data  points. 


APPENDIX  D 

ANALYTIC  REPRESENTATION  OF  STEADY  STATE 
ONE  SPECIES  TWO  BODY  DECAY 


We  present  here  the  derivation  of  an  analytic  representation  of  the 
steady  state  solution  of  one  species  two  body  decay.  The  asymptotic 
limits  we  seek  to  recover  are  three:  (1)  pure  diffusion  as  k/S  •>  0,  (2) 
pure  localized  kinetics  as  k/3  » within  the  central  core,  and  (3) 

source  independent  diffusion  kinetic  balance  in  the  limb. 

The  method  which  suggests  itself  is  to  gate  between  core  (pure 
localized  kinetics)  and  limb  (source  independent  diffusion  kinetic 
balance)  for  k/ 3 y 0,  and  then  to  gate  this  expression  over  to  the  pure 
diffusion  limit  as  k/ 3 0.  Since  the  core  is  self-gating  in  k/3,  it 

is  required  to  gate  the  limb  from  the  1/r^  diffusion  kinetic  limit  over 
to  the  1/r  pure  diffusion  limit. 

In  actual  fact  it  proved  convenient  to  go  about  this  latter  step 
backwards:  one  gates  to  the  pure  diffusion  limit,  but  not  just  the  pure 

diffusion  limb  limit,  but  the  global  pure  diffusion  limit,  and  then 
corrects  for  it  in  the  limb  for  k/3  > 0,  via  a gate  in  r. 

The  form  is  given  by 


(D.l) 


where 


(D.2) 
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'’kin 


where 


Y = 


and 


gS/2< 

^ 3^.2  J-B 


= B/(3^  + K^). 


(D.3) 


(D.4) 


(D.5) 


The  two  choices  for  the  gates  are 


(1)  f - g - 

e + K 0 + K 


and 


(2)  f = 


}2p 


,2p 


S^P  * K^P  ’ ® ■ j2p  + ,2p 


And,  finally. 


I +^[(r‘'  + r I)'/*!  - r 1 

g J O o 

r = r'/B. 
o o 


(D.6) 


(D.7) 


(D.8) 


where 
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Certain  features  of  this  form  deserve  explanation. 

1.  In  the  limit  k = 0,  G = For  k > 0,  G - 2B/K^r  as 

r oo,  taking  advantage  of  the  constant  term,  2,  within 

brackets  in  G,.c. 

dit 

2.  Since  the  global  diffusion  solution  is  used,  which  is  slope 
regularized  at  r=0,  and  since  the  factor  F is  slope  regularized 

(for  q > 1),  their  product  is  also.  By  construction,  the 
difference  of  exponentials  is  slope  regularized  at  r=0, 
irrespective  of  the  value  of  y,  wherefore  the  total  sum  is  also. 

3.  Allowing  the  opposite  gates  in  ic/B,  f and  g,  to  remain 

completely  general  for  the  moment,  y is  constructed  to  regain 
the  result,  C(0)  = B/(3  + k^).  The  parameter  y does  however 

have  a further  constraint  on  it,  namely  that  y > 3/2.  This  is 
met  by  setting  p > 1 in  f and  g for  either  choice  of  gate. 


APPENDIX  E 

DIFFUSION  KINETICS  FOR  v BODY  DECAY 
As  an  aside  of  theoretical  interest  we  consider  diffusion  kinetics 
for  v-body  decay  of  one  species  for  a continuous  source  term,  turned  on 
at  t = 0.  The  equation  is  given  by 


8c 

at 


D a , 2 ac.  , V 

r 


(E.l) 


For  the  spatial  dependence  of  the  source  term,  we  use  the  now  familiar 
exponential  form,  y(r)  = Ae"^’'  = ( B^Z/8Tr)e"^’^,  where  Z equals  the 
spatially  integrated  source.  In  the  limit  B y(r)  > (Z/4Tr)6(r). 

The  more  general  form,  w(r)  = ( b'^^Z/AttPC \H-3) )r''e~^’^,  where  v > 0, 
also  approaches  (Z/A-n)S(r)  as  B ->■  «. 

The  problem  in  its  full  time  dependent  context,  specifying  zero 
initial  conditions,  is  well  posed.  The  solution,  c(r,t),  at  earliest 
times  is  given  by 


c(r,t)  ~ y(r)t  (E.2) 

Although  the  existence,  uniqueness  and  positive  definiteness  of  the 
solution  at  all  times,  t > 0,  may  be  proved  by  rigorous  argument,  they 
are  clear  from  a physical  understanding  of  the  equation. 

We  consider  in  particular  the  steady  state  limit.  Setting  c = g/r, 
Eq.  (E.l)  becomes 
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g"  (f)'’r  = -(3^Z/8TTD)re  (e.3) 

Since  the  steady  state  limit  is  the  closure  of  the  time  dependent 
solution,  one  may  only  argue  that  it  be  nonnegative  definite.  Consider 
now  either  the  long  range  limit  of  the  solution  for  given  B or  the 
solution  at  all  r > 0 in  the  limit  B Equation  (E.3)  then  becomes 


r = 0 


(E.4) 


for  r > 0.  Substituting  g = as  a trial  solution,  one  obtains  the 

requirements 


a = 


v-3 

v-1 


and 


®o  = 

(v-l)^(vk/D) 

« 


(E.5) 


These  obviously  do  not  apply  to  the  linear  cases  v = 0 and  v = 1 
which  are  covered  separateiy  by  the  well  known  solutions,  c(r)  = Z/Dr  in 
the  first  case  and  c(r)  = Ze  '^^/Dr,  where  k = /k/D,  in  the  second. 

Note  that  c(r)  in  both  of  these  cases  is  homogeneous  in  the  source 
strength,  Z. 

1^  ^11  the  other  cases , the  noniinear  ones , this  long  range 
solution  (or  the  solution  for  r > 0 in  the  limit  B ^ is  independent 
of  the  strength  of  the  source,  and  is  uniquely  determined  by  the 
equation  without  boundary  conditions,  with  the  exception  that  we  require 
the  solution  to  decrease  to  zero  as  r -»■  ».  Now,  looking  at  these 
requirements,  we  see  that,  given  the  trial  solution  form,  the  nonzero 
solution  is  positive  definite  only  for  0 < v < 3,  degenerates  to  the 
zero  solution  for  v = 3 and  is  in  general  complex  for  v > 3.  We 
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note,  however,  that  for  0 < v < 1,  the  solution  does  not  decrease  to 
zero  as  r + 0.  This  is  unacceptable.  We  note,  curiously  enough,  for 
the  complex  cases  (1.  e.  v > 3)  that  the  modulus  of  the  solution  does 
decrease  to  zero  as  r -»•  0. 

Arguing  from  the  analyticity  of  the  form  of  the  equation,  we 
conclude  that  the  only  physically  acceptable  solution  which  exists  in 
the  intervals  0 < v < 1 and  3 < v is  the  zero  solution.  Why  this  is  so 
appears  as  something  of  a mystery  from  the  isolated  vantage  point  of  the 
steady  state  problem  by  itself.  It  should  become  clear  upon  considering 
the  full  time  dependent  solution.  In  this  regard  the  pure  kinetic 
results  for  general  v may  prove  relevant.  It  would  be  an  interesting 
direction  in  which  to  take  this  thesis  work — namely,  to  expand  the 
capabilities  of  the  computer  code,  DIFKIN,  to  calculate  for  general 
v-body  kinetic  terms. 

It  is  possible  to  generalize  the  inverse  length,  k,  defined  in 
the  two-body  case  to  the  general  v-body  case  (v  * 0,  ^ 1) . In  the 
pure  kinetic  limit,  c(0)  = ^A/vk.  If  we  there  set  c(0)  = B/k^  where 
B = A/D  as  was  done  in  the  two-body  case,  then  we  define  an  inverse 
length  associated  with  the  problem. 


.DA  .l/v.1/2 


(E.6) 


It  may  prove  fruitful  to  consider  the  envelope  of  solutions  for  constant 
K for  fixed  A and  D.  It  is  conjectured  that  the  form  B/(B^  + k^) 
should  again  be  a reasonably  good  approximation  to  c(0).  This 
approximate  form  looks  like  a first  order  result  in  a Fourier  transform 
space.  It  should  prove  interesting  to  pursue  this  work  further. 


APPENDIX  F 
COMPUTER  WORK 

In  this  appendix  a library  of  the  major  computer  programs  and  data 
files  written  or  adapted  or  compiled  for  or  implemented  in  this  thesis 
are  given.  They  are  listed  in  the  order  of  the  chapter  or  appendix  they 
refer  to,  and  are  given  in  the  format  of  TCP  files,  in  many  cases  job 
files  with  /*include's  and  JCL.  The  meaning  is  generally  clear  from 
context  and,  in  most  cases,  clearer  by  virtue  of  displaying  in  a rather 
convenient  manner  program  and  job  organization. 

Each  chapter  section  is  so  indicated  by  a header.  Each  TCP  file  is 
preceded  by  a header,  giving  its  name.  This  is  important,  since  this  is 
how  it  is  referenced  in  any  /*include's  appearing  in  any  job  files  which 
may  include  it.  Adequate  documentation  is  generally  given  in  each 
program,  explaining  the  purpose  of  the  program  and  what  is  being  done  in 
it. 
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REAO( 2,207)  NJ (I) , ISPECFI I ) , ( NJ J ( I , NN ) ,NN=1 , 19 ) 
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